OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
alew7.F File Reference
#include "implicit_f.inc"
#include "task_c.inc"
#include "comlock.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine alew7 (x, v, w, ms, nale, nodft, nodlt, weight, numnod, dt1, sx, sv, sw, nspmd)

Function/Subroutine Documentation

◆ alew7()

subroutine alew7 ( dimension(3,sx/3), intent(inout) x,
dimension(3,sv/3), intent(in) v,
dimension(3,sw/3), intent(inout) w,
dimension(numnod), intent(in) ms,
integer, dimension(numnod), intent(in) nale,
integer, intent(in) nodft,
integer, intent(in) nodlt,
integer, dimension(numnod), intent(in) weight,
integer, intent(in) numnod,
intent(in) dt1,
integer, intent(in) sx,
integer, intent(in) sv,
integer, intent(in) sw,
integer, intent(in) nspmd )

Definition at line 42 of file alew7.F.

46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE ale_mod
50 USE spmd_exch_flow_tracking_data_mod
51 USE spmd_exch_flow_tracking_data2_mod
52 USE spmd_exch_flow_tracking_data3_mod
53 USE spmd_exch_flow_tracking_data4_mod
54C-----------------------------------------------
55C D e s c r i p t i o n
56C-----------------------------------------------
57C Compute Grid velocities for /ALE/GRID/FLOW-TRACKING
58C
59C INPUT/OUTPUT
60C X,D,V are allocated to SX,SD,DV=3*(NUMNOD_L+NUMVVOIS_L)
61C in grid subroutine it may needed to access nodes which
62C are connected to a remote elem. They are sored in X(1:3,NUMNOD+1:)
63C Consequently X is defined here X(3,SX/3) instead of X(3,NUMNOD) as usually
64C WEIGHT is a tag (0 or 1) which ensure a unique contribution of a given nodal value when sum is calculated
65C over all SPMD domain. (otherwise nodes at common boundary with another domain will have multiple contributions)
66C
67C A AVERAGED VALUE IS COMPUTED IN THE MASS FLOW : VMEAN(1:2)
68C THIS VALUE IS APPLIED TO THE GRID VELOCITY
69C SUM_MOM(1) : is Sum(m[i]*vx[i], i=1..numnod)
70C SUM_MOM(2) : is Sum(m[i]*vy[i], i=1..numnod)
71C SUM_MOM(3) : is Sum(m[i]*vz[i], i=1..numnod)
72C SUM_MOM(4) : is Sum(m[i] , i=1..numnod)
73C
74C VMEAN(1) = SUM_MOM(1) / SUM_MOM(4)
75C VMEAN(2) = SUM_MOM(1) / SUM_MOM(4)
76C VMEAN(3) = SUM_MOM(1) / SUM_MOM(4)
77C
78C SMP :
79C SUM_MOM(1:4) are the cumulative from i=NODFT,NODLT
80C SUM_MOM_L(1:4) is then the cumulative result from i=1,NUMNOD
81C
82C SPMD :
83C SUM_MOM_L(1:4) is calculated on each domain
84C then exchanged
85C finally each domain can deduce the complete cumulative result
86C
87C PARITH/ON :
88C subroutine FOAT_TO_6_FLOAT is used for this purpose
89C
90C-----------------------------------------------
91C I m p l i c i t T y p e s
92C-----------------------------------------------
93#include "implicit_f.inc"
94C-----------------------------------------------
95C C o m m o n B l o c k s
96C-----------------------------------------------
97#include "task_c.inc"
98#include "comlock.inc"
99C-----------------------------------------------
100C D u m m y A r g u m e n t s
101C-----------------------------------------------
102! SPMD CASE : SX >= 3*NUMNOD (SX = 3*(NUMNOD_L+NRCVVOIS_L))
103! X(1:3,1:NUMNOD) : local nodes
104! (1:3, NUMNOD+1:) additional nodes (also on adjacent domains but connected to the boundary of the current domain)
105! idem with D(SD), and V(SV)
106C-----------------------------------------------
107 INTEGER, INTENT(IN) :: NSPMD,NUMNOD,SX,SV,SW
108 INTEGER, INTENT(IN) :: NALE(NUMNOD), NODFT, NODLT
109 INTEGER, INTENT(IN) :: WEIGHT(NUMNOD)
110 my_real, INTENT(INout) :: x(3,sx/3)
111 my_real, INTENT(IN) :: v(3,sv/3), ms(numnod)
112 my_real, INTENT(INOUT) :: w(3,sw/3)
113 my_real, INTENT(IN) :: dt1
114C-----------------------------------------------
115C L o c a l V a r i a b l e s
116C-----------------------------------------------
117 INTEGER :: I,JJ
118 my_real :: sum_ms,ms_node, ksi
119 my_real :: sum_mom(3),sum_cog(3) !1:m*vx, 2:m*vy, 3:m*vz
120 my_real :: sum_m
121 my_real :: sum_itm(6)
122 my_real :: vmean(3),ld(6),lw(3),cog(3),itm(6),ld_b(3,3),ld_bp(3,3),ld_bp_b(3,3),ld_tot_b(3,3)
123 my_real :: p_b_bp(3,3)
124 my_real :: eigenval(3),eigenvec(3,3)
125 my_real :: x_min_max(6),x_min_max_grid(6)
126 my_real :: r11,r12,r13, r21,r22,r23, r31,r32,r33 ! transition matrix
127 my_real :: x_trans(3) ! vector after transition matrix
128 my_real :: xx,yy,zz
129 my_real :: ld_norm
130 my_real :: scale_def,scale_rot
131 my_real :: ratio(3), dw(3)
132 my_real :: beta(6),beta0(6)
133 my_real :: fac(3)
134 my_real :: ms_elem_mean
135 LOGICAL :: lTENSOR, lDEF,lROT
136 LOGICAL :: HAS_ALE_NODE, HAS_FLOW_NODE
137C-----------------------------------------------
138C S o u r c e L i n e s
139C-----------------------------------------------
140 cog(1:3) = zero
141 !1.---user parameters
142 !-----------------------------------
143 ltensor = .false.
144 ldef=.false.
145 lrot=.false.
146 scale_def = ale%GRID%ALPHA !user parameter
147 scale_rot = ale%GRID%GAMMA !user parameter
148 IF(int(ale%GRID%VGX) == 1)ldef=.true. !enable mesh deformation
149 IF(int(ale%GRID%VGY) == 1)lrot=.true. !enable mesh rotation
150 IF(ldef .OR. lrot)ltensor=.true.!request for mean tensor calculation
151
152 !2.---AVERAGED DEFORMATION TENSOR
153 !-----------------------------------
154 !mean value is computed here
155 ! cumulative mass value done in sfor3.F for 3d solid and qforc2.F for quads
156 ! ALE%GRID%flow_tracking_data%.. is a shared data structure for SMP (omp single or lock on/off)
157!$OMP SINGLE
158 IF(ltensor)THEN
159 IF(nspmd > 1)THEN
160 CALL spmd_exch_flow_tracking_data(ale%GRID%flow_tracking_data, nspmd)
161 ENDIF
162 ale%GRID%flow_tracking_data%EP(1:9) = ale%GRID%flow_tracking_data%EP(1:9) / ale%GRID%flow_tracking_data%SUM_M
163 !STRAIN RATE TENSOR
164 ale%GRID%flow_tracking_data%LD(1) = ale%GRID%flow_tracking_data%EP(1) !xx
165 ale%GRID%flow_tracking_data%LD(2) = ale%GRID%flow_tracking_data%EP(2) !yy
166 ale%GRID%flow_tracking_data%LD(3) = ale%GRID%flow_tracking_data%EP(3) !zz
167 ale%GRID%flow_tracking_data%LD(4) = half*(ale%GRID%flow_tracking_data%EP(4)+ale%GRID%flow_tracking_data%EP(7)) !xy
168 ale%GRID%flow_tracking_data%LD(5) = half*(ale%GRID%flow_tracking_data%EP(5)+ale%GRID%flow_tracking_data%EP(8)) !xz
169 ale%GRID%flow_tracking_data%LD(6) = half*(ale%GRID%flow_tracking_data%EP(6)+ale%GRID%flow_tracking_data%EP(9)) !yz
170 !SPIN TENSOR
171 ale%GRID%flow_tracking_data%LW(1) = half*(ale%GRID%flow_tracking_data%EP(4)-ale%GRID%flow_tracking_data%EP(7)) !z
172 ale%GRID%flow_tracking_data%LW(2) = half*(ale%GRID%flow_tracking_data%EP(5)-ale%GRID%flow_tracking_data%EP(8)) !y
173 ale%GRID%flow_tracking_data%LW(3) = half*(ale%GRID%flow_tracking_data%EP(6)-ale%GRID%flow_tracking_data%EP(9)) !x
174 ENDIF
175!$OMP END SINGLE
176
177 CALL my_barrier
178
179 !4.---MEAN VELOCITY & CENTER OF MASS
180 !-----------------------------------
181 ale%GRID%flow_tracking_data%MOM_L(1:3) = zero
182 ale%GRID%flow_tracking_data%COG_L(1:3) = zero
183 ale%GRID%flow_tracking_data%SUM_M = zero !used for elem mass in sforc3, used below for nodal mass
184 sum_mom(1:3) = zero
185 sum_cog(1:3) = zero
186 sum_m = zero
187 DO i = nodft, nodlt
188 IF(iabs(nale(i)) == 1)THEN
189 ms_node = ms(i)*weight(i) !WEIGHT(1:NUMNOD) ENSURE A SINGLE CONTRIBUTION FROM ALL DOMAIN
190 sum_mom(1) = sum_mom(1) + ms_node*v(1,i)
191 sum_mom(2) = sum_mom(2) + ms_node*v(2,i)
192 sum_mom(3) = sum_mom(3) + ms_node*v(3,i)
193 sum_m = sum_m + ms_node
194 sum_cog(1) = sum_cog(1) + ms_node*x(1,i)
195 sum_cog(2) = sum_cog(2) + ms_node*x(2,i)
196 sum_cog(3) = sum_cog(3) + ms_node*x(3,i)
197 ENDIF
198 ENDDO
199 CALL my_barrier
200#include "lockon.inc"
201 !assembly of global values for SMP case
202 ale%GRID%flow_tracking_data%MOM_L(1) = ale%GRID%flow_tracking_data%MOM_L(1) + sum_mom(1)
203 ale%GRID%flow_tracking_data%MOM_L(2) = ale%GRID%flow_tracking_data%MOM_L(2) + sum_mom(2)
204 ale%GRID%flow_tracking_data%MOM_L(3) = ale%GRID%flow_tracking_data%MOM_L(3) + sum_mom(3)
205 ale%GRID%flow_tracking_data%SUM_M = ale%GRID%flow_tracking_data%SUM_M + sum_m
206 ale%GRID%flow_tracking_data%COG_L(1) = ale%GRID%flow_tracking_data%COG_L(1) + sum_cog(1)
207 ale%GRID%flow_tracking_data%COG_L(2) = ale%GRID%flow_tracking_data%COG_L(2) + sum_cog(2)
208 ale%GRID%flow_tracking_data%COG_L(3) = ale%GRID%flow_tracking_data%COG_L(3) + sum_cog(3)
209#include "lockoff.inc"
210 CALL my_barrier
211 !SPMD EXCHANGE 'MOM' and 'COG'
212
213!$OMP SINGLE
214 IF(nspmd > 1)THEN
215 CALL spmd_exch_flow_tracking_data2(ale%GRID%flow_tracking_data, nspmd)
216 ENDIF
217!$OMP END SINGLE
218
219 vmean(1:3) = zero
220 sum_ms = ale%GRID%flow_tracking_data%SUM_M
221 IF(sum_ms > em20)THEN
222 vmean(1) = ale%GRID%flow_tracking_data%MOM_L(1)/sum_ms
223 vmean(2) = ale%GRID%flow_tracking_data%MOM_L(2)/sum_ms
224 vmean(3) = ale%GRID%flow_tracking_data%MOM_L(3)/sum_ms
225 cog(1) = ale%GRID%flow_tracking_data%COG_L(1)/sum_ms
226 cog(2) = ale%GRID%flow_tracking_data%COG_L(2)/sum_ms
227 cog(3) = ale%GRID%flow_tracking_data%COG_L(3)/sum_ms
228 END IF
229
230
231 !5.---INERTIA TENSOR MATRIX (ITM)
232 !-----------------------------------
233 ale%GRID%flow_tracking_data%ITM_L(1:6) = zero
234 CALL my_barrier
235 sum_itm(1:6) = zero
236 DO i = nodft, nodlt
237 IF(iabs(nale(i)) == 1)THEN
238 xx = x(1,i)-cog(1)
239 yy = x(2,i)-cog(2)
240 zz = x(3,i)-cog(3)
241 ms_node = ms(i)*weight(i) !WEIGHT(1:NUMNOD) ENSURE A SINGLE CONTRIBUTION FROM ALL DOMAIN
242 sum_itm(1) = sum_itm(1) + ms_node*(yy*yy+zz*zz)
243 sum_itm(2) = sum_itm(2) + ms_node*(xx*xx+zz*zz)
244 sum_itm(3) = sum_itm(3) + ms_node*(xx*xx+yy*yy)
245 sum_itm(4) = sum_itm(4) - ms_node*(xx*yy)
246 sum_itm(5) = sum_itm(5) - ms_node*(yy*zz)
247 sum_itm(6) = sum_itm(6) - ms_node*(xx*zz)
248 ENDIF
249 ENDDO
250#include "lockon.inc"
251 !assembly of global values for SMP case
252 ale%GRID%flow_tracking_data%ITM_L(1) = ale%GRID%flow_tracking_data%ITM_L(1) + sum_itm(1)
253 ale%GRID%flow_tracking_data%ITM_L(2) = ale%GRID%flow_tracking_data%ITM_L(2) + sum_itm(2)
254 ale%GRID%flow_tracking_data%ITM_L(3) = ale%GRID%flow_tracking_data%ITM_L(3) + sum_itm(3)
255 ale%GRID%flow_tracking_data%ITM_L(4) = ale%GRID%flow_tracking_data%ITM_L(4) + sum_itm(4)
256 ale%GRID%flow_tracking_data%ITM_L(5) = ale%GRID%flow_tracking_data%ITM_L(5) + sum_itm(5)
257 ale%GRID%flow_tracking_data%ITM_L(6) = ale%GRID%flow_tracking_data%ITM_L(6) + sum_itm(6)
258#include "lockoff.inc"
259 CALL my_barrier
260 !SPMD EXCHANGE 'ITM'
261!$OMP SINGLE
262 IF(nspmd > 1)THEN
263 CALL spmd_exch_flow_tracking_data3(ale%GRID%flow_tracking_data, nspmd)
264 ENDIF
265!$OMP END SINGLE
266 sum_ms = ale%GRID%flow_tracking_data%SUM_M
267 IF(sum_ms > em20)THEN
268 itm(1) = ale%GRID%flow_tracking_data%ITM_L(1)/sum_ms
269 itm(2) = ale%GRID%flow_tracking_data%ITM_L(2)/sum_ms
270 itm(3) = ale%GRID%flow_tracking_data%ITM_L(3)/sum_ms
271 itm(4) = ale%GRID%flow_tracking_data%ITM_L(4)/sum_ms
272 itm(5) = ale%GRID%flow_tracking_data%ITM_L(5)/sum_ms
273 itm(6) = ale%GRID%flow_tracking_data%ITM_L(6)/sum_ms
274 END IF
275
276
277 !6.---EIGENVECTORS - INERTIA TENSOR
278 !-----------------------------------
279 CALL valpvec(itm,eigenval,eigenvec,1)
280 !-----------------------------------
281 IF(dt1 == zero)THEN
282 ale%GRID%flow_tracking_data%EIGENVEC(1:3,1)=eigenvec(1:3,1)
283 ale%GRID%flow_tracking_data%EIGENVEC(1:3,2)=eigenvec(1:3,2)
284 ale%GRID%flow_tracking_data%EIGENVEC(1:3,3)=eigenvec(1:3,3)
285 ENDIF
286 eigenvec(1:3,1) = ale%GRID%flow_tracking_data%EIGENVEC(1:3,1)
287 eigenvec(1:3,2) = ale%GRID%flow_tracking_data%EIGENVEC(1:3,2)
288 eigenvec(1:3,3) = ale%GRID%flow_tracking_data%EIGENVEC(1:3,3)
289
290
291 !7.---TRANSITION MATRIX (transposed)
292 !-----------------------------------
293 r11=ale%GRID%flow_tracking_data%EIGENVEC(1,1)
294 r21=ale%GRID%flow_tracking_data%EIGENVEC(1,2)
295 r31=ale%GRID%flow_tracking_data%EIGENVEC(1,3)
296
297 r12=ale%GRID%flow_tracking_data%EIGENVEC(2,1)
298 r22=ale%GRID%flow_tracking_data%EIGENVEC(2,2)
299 r32=ale%GRID%flow_tracking_data%EIGENVEC(2,3)
300
301 r13=ale%GRID%flow_tracking_data%EIGENVEC(3,1)
302 r23=ale%GRID%flow_tracking_data%EIGENVEC(3,2)
303 r33=ale%GRID%flow_tracking_data%EIGENVEC(3,3)
304
305
306 !8.---MEAN MASS DENSITY
307 !-----------------------------------
308 CALL my_barrier
309
310 IF(dt1 == zero)THEN
311 ms_elem_mean = ep20
312 IF(ale%GRID%flow_tracking_data%NUM_ELEM_ALE > 0)THEN
313 ms_elem_mean = ale%GRID%flow_tracking_data%SUM_M / (one*ale%GRID%flow_tracking_data%NUM_ELEM_ALE)
314 ENDIF
315 ale%GRID%flow_tracking_data%MS_ELEM_MEAN_0 = ms_elem_mean
316 ENDIF
317 ms_elem_mean = ale%GRID%flow_tracking_data%MS_ELEM_MEAN_0
318
319
320 !9.---ENCOMPASSING BOX (FLOW & WHOLE GRID)
321 !-----------------------------------
322 ale%GRID%flow_tracking_data%X_MIN_MAX(1:3) = ep20 !FLOW min 1,2,3
323 ale%GRID%flow_tracking_data%X_MIN_MAX(4:6) = -ep20 !FLOW max 1,2,3
324 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(1:3) = ep20 ! FULL GRID min 1,2,3
325 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(4:6) = -ep20 ! FULL GRID max 1,2,3
326 CALL my_barrier
327
328 x_min_max(1:3) = ep20
329 x_min_max(4:6) = -ep20
330 x_min_max_grid(1:3) = ep20
331 x_min_max_grid(4:6) = -ep20
332
333 has_ale_node = .false.
334 has_flow_node = .false.
335
336 DO i = nodft, nodlt
337 IF(iabs(nale(i)) == 1)THEN
338 has_ale_node = .true.
339 ! apply transition matrix (local basis)
340 x_trans(1) = r11*(x(1,i)-cog(1))+r12*(x(2,i)-cog(2))+r13*(x(3,i)-cog(3))
341 x_trans(2) = r21*(x(1,i)-cog(1))+r22*(x(2,i)-cog(2))+r23*(x(3,i)-cog(3))
342 x_trans(3) = r31*(x(1,i)-cog(1))+r32*(x(2,i)-cog(2))+r33*(x(3,i)-cog(3))
343 !FULL GRID encompassing box (local basis)
344 x_min_max_grid(1) = min(x_min_max_grid(1),x_trans(1))
345 x_min_max_grid(2) = min(x_min_max_grid(2),x_trans(2))
346 x_min_max_grid(3) = min(x_min_max_grid(3),x_trans(3))
347 x_min_max_grid(4) = max(x_min_max_grid(4),x_trans(1))
348 x_min_max_grid(5) = max(x_min_max_grid(5),x_trans(2))
349 x_min_max_grid(6) = max(x_min_max_grid(6),x_trans(3))
350 !MAIN FLOW encompassing box(local basis)
351 ms_node = ms(i)*weight(i)
352 ratio(1) = ms_node / ms_elem_mean !MS_ELEM_MEAN_0 ! with ratio : focus only on main flow ! w/o ratio : all coordinates
353 IF(ratio(1) > one)THEN
354 has_flow_node = .true.
355 x_min_max(1) = min(x_min_max(1),x_trans(1))
356 x_min_max(2) = min(x_min_max(2),x_trans(2))
357 x_min_max(3) = min(x_min_max(3),x_trans(3))
358 x_min_max(4) = max(x_min_max(4),x_trans(1))
359 x_min_max(5) = max(x_min_max(5),x_trans(2))
360 x_min_max(6) = max(x_min_max(6),x_trans(3))
361 ENDIF
362 ENDIF
363 ENDDO
364 !gathering values for SMP
365#include "lockon.inc"
366 !assembly of global values for SMP case
367 !---FLOW
368 IF(has_flow_node)THEN
369 ale%GRID%flow_tracking_data%X_MIN_MAX(1) = min(ale%GRID%flow_tracking_data%X_MIN_MAX(1), x_min_max(1))
370 ale%GRID%flow_tracking_data%X_MIN_MAX(2) = min(ale%GRID%flow_tracking_data%X_MIN_MAX(2), x_min_max(2))
371 ale%GRID%flow_tracking_data%X_MIN_MAX(3) = min(ale%GRID%flow_tracking_data%X_MIN_MAX(3), x_min_max(3))
372 ale%GRID%flow_tracking_data%X_MIN_MAX(4) = max(ale%GRID%flow_tracking_data%X_MIN_MAX(4), x_min_max(4))
373 ale%GRID%flow_tracking_data%X_MIN_MAX(5) = max(ale%GRID%flow_tracking_data%X_MIN_MAX(5), x_min_max(5))
374 ale%GRID%flow_tracking_data%X_MIN_MAX(6) = max(ale%GRID%flow_tracking_data%X_MIN_MAX(6), x_min_max(6))
375 ENDIF
376 !---FULL GRID
377 IF(has_ale_node)THEN
378 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(1) = min(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(1), x_min_max_grid(1))
379 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(2) = min(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(2), x_min_max_grid(2))
380 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(3) = min(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(3), x_min_max_grid(3))
381 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(4) = max(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(4), x_min_max_grid(4))
382 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(5) = max(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(5), x_min_max_grid(5))
383 ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(6) = max(ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(6), x_min_max_grid(6))
384 ENDIF
385#include "lockoff.inc"
386 CALL my_barrier
387 !gathering values for SPMD
388!$OMP SINGLE
389 IF(nspmd > 1)THEN
390 CALL spmd_exch_flow_tracking_data4(ale%GRID%flow_tracking_data, nspmd)
391 ENDIF
392!$OMP END SINGLE
393 CALL my_barrier
394
395
396 !10. RATIO (to detect border proximity)
397 !-----------------------------------
398 !local axis 1 - BETA_Xmin,BETA_Xmax
399 beta(1)=ale%GRID%flow_tracking_data%X_MIN_MAX(1)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(1)
400 beta(4)=ale%GRID%flow_tracking_data%X_MIN_MAX(4)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(4)
401 !local axis 2 - BETA_Ymin,BETA_Ymax
402 beta(2)=ale%GRID%flow_tracking_data%X_MIN_MAX(2)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(2)
403 beta(5)=ale%GRID%flow_tracking_data%X_MIN_MAX(5)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(5)
404 !local axis 3 - BETA_Zmin,BETA_Zmax
405 beta(3)=ale%GRID%flow_tracking_data%X_MIN_MAX(3)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(3)
406 beta(6)=ale%GRID%flow_tracking_data%X_MIN_MAX(6)/ale%GRID%flow_tracking_data%X_MIN_MAX_GRID(6)
407 ! store initial value
408 IF(dt1 == zero)THEN
409 ale%GRID%flow_tracking_data%BETA0(1) = beta(1)
410 ale%GRID%flow_tracking_data%BETA0(2) = beta(2)
411 ale%GRID%flow_tracking_data%BETA0(3) = beta(3)
412 ale%GRID%flow_tracking_data%BETA0(4) = beta(4)
413 ale%GRID%flow_tracking_data%BETA0(5) = beta(5)
414 ale%GRID%flow_tracking_data%BETA0(6) = beta(6)
415 END IF
416 CALL my_barrier
417 !retrieve initial ratio
418 beta0(1) = ale%GRID%flow_tracking_data%BETA0(1)
419 beta0(2) = ale%GRID%flow_tracking_data%BETA0(2)
420 beta0(3) = ale%GRID%flow_tracking_data%BETA0(3)
421 beta0(4) = ale%GRID%flow_tracking_data%BETA0(4)
422 beta0(5) = ale%GRID%flow_tracking_data%BETA0(5)
423 beta0(6) = ale%GRID%flow_tracking_data%BETA0(6)
424
425
426 !11.---NORM OF STRAIN RATE TENSOR (global basis)
427 !-----------------------------------
428 CALL my_barrier
429 !local working array
430 ld(1:6)=ale%GRID%flow_tracking_data%LD(1:6)
431 lw(1:3)=ale%GRID%flow_tracking_data%LW(1:3)
432 ale%GRID%flow_tracking_data%SUM_M = zero
433
434 ld_b(1,1:3) = (/ld(1),ld(4),ld(5)/)
435 ld_b(2,1:3) = (/ld(4),ld(2),ld(6)/)
436 ld_b(3,1:3) = (/ld(5),ld(6),ld(3)/)
437
438 ld_norm = zero
439 ld_norm = ld_norm + abs(ld_b(1,1)) + abs(ld_b(1,2)) + abs(ld_b(1,3))
440 ld_norm = ld_norm + abs(ld_b(2,1)) + abs(ld_b(2,2)) + abs(ld_b(2,3))
441 ld_norm = ld_norm + abs(ld_b(3,1)) + abs(ld_b(3,2)) + abs(ld_b(3,3))
442 ld_norm = ld_norm / nine
443
444 CALL my_barrier
445 IF( ale%GRID%flow_tracking_data%LD_NORM < ld_norm)THEN
446 ale%GRID%flow_tracking_data%LD_NORM = ld_norm
447 ENDIF
448 ld_norm = ale%GRID%flow_tracking_data%LD_NORM
449
450
451 !12.---CORRECTION TERM : LD_TOTAL = LD_AVERAGE + LD_CORRECTION
452 !-----------------------------------
453 ! B : global basis
454 ! B' : local basis (Inertia Tensor Matrix)
455 !transition matrix P_[B->B']
456 p_b_bp(1:3,1)=eigenvec(1:3,1)
457 p_b_bp(1:3,2)=eigenvec(1:3,2)
458 p_b_bp(1:3,3)=eigenvec(1:3,3)
459
460 !factor margin in basis B'
461 ! fac > 1 means that border is too near the flow
462 fac(1) = max(one, beta(1)/beta0(1), beta(4)/beta0(4))
463 fac(2) = max(one, beta(2)/beta0(2), beta(5)/beta0(5))
464 fac(3) = max(one, beta(3)/beta0(3), beta(6)/beta0(6))
465
466 ! amplification with zero derivative at 1.000
467 ksi = ep02
468 fac(1) = half*ld_norm *(one + ksi*(fac(1)-one)**2)
469 fac(2) = half*ld_norm *(one + ksi*(fac(2)-one)**2)
470 fac(3) = half*ld_norm *(one + ksi*(fac(3)-one)**2)
471
472 ! correction term (strain rate tensor) in local basis B'
473 ld_bp(1,1:3) = (/ fac(1) , zero , zero /)
474 ld_bp(2,1:3) = (/ zero , fac(2) , zero /)
475 ld_bp(3,1:3) = (/ zero , zero , fac(3) /)
476
477 ! correction term (strain rate tensor) in in global basis B
478 ld_bp_b = matmul( p_b_bp, matmul(ld_bp, transpose(p_b_bp)) )
479
480
481 !13.---TOTAL STRAIN RATE TENSOR
482 !-----------------------------------
483 !LD_B : averaged strain rate tensor in global basis
484 !LD_Bp_B : correction term computed in local basis and then written in global basis
485 !LD_TOT_B : total strain rate tensor in global basis
486 ld_tot_b = ld_b + ld_bp_b
487 !working array
488 ld(1)=ld_tot_b(1,1)
489 ld(2)=ld_tot_b(2,2)
490 ld(3)=ld_tot_b(3,3)
491 ld(4)=ld_tot_b(1,2)
492 ld(5)=ld_tot_b(1,3)
493 ld(6)=ld_tot_b(2,3)
494
495
496 !14.---UPDATE GRID
497 !-----------------------------------
498 DO i = nodft, nodlt
499 IF(iabs(nale(i)) == 1) THEN
500 ! GRID TRANSLATION
501 w(1,i)=vmean(1)
502 w(2,i)=vmean(2)
503 w(3,i)=vmean(3)
504
505 ! coordinates relative to center of mass
506 xx = (x(1,i)-cog(1))
507 yy = (x(2,i)-cog(2))
508 zz = (x(3,i)-cog(3))
509
510 !GRID DEFORMATION
511 IF(ldef)THEN
512 !increment in global basis (from strain rate tensor + correction term)
513 dw(1) = scale_def*(ld(1)*xx+ld(4)*yy+ld(5)*zz)
514 dw(2) = scale_def*(ld(4)*xx+ld(2)*yy+ld(6)*zz)
515 dw(3) = scale_def*(ld(5)*xx+ld(6)*yy+ld(3)*zz)
516 ! update grid velocities
517 w(1,i) =w(1,i)+ dw(1)
518 w(2,i) =w(2,i)+ dw(2)
519 w(3,i) =w(3,i)+ dw(3)
520 ENDIF
521
522 !GRID ROTATION
523 IF(lrot)THEN
524 w(1,i) = w(1,i) + scale_rot*(+lw(1)*yy+lw(2)*zz)
525 w(2,i) = w(2,i) + scale_rot*(-lw(1)*xx+lw(3)*zz)
526 w(3,i) = w(3,i) + scale_rot*(-lw(2)*xx-lw(3)*yy)
527 ENDIF
528
529 ELSEIF(nale(i) == 0)THEN
530 ! lagrangian framework
531 w(1,i)=v(1,i)
532 w(2,i)=v(2,i)
533 w(3,i)=v(3,i)
534 ELSE
535 ! eulerian framework
536 w(1,i)=zero
537 w(2,i)=zero
538 w(3,i)=zero
539 ENDIF
540 ENDDO
541
542
543 !14.---LOCAL AXIS UPDATE (rotation)
544 !-----------------------------------
545 IF(lrot)THEN
546 !rotate local basis
547 DO jj=1,3
548 xx = eigenvec(1,jj)
549 yy = eigenvec(2,jj)
550 zz = eigenvec(3,jj)
551 eigenvec(1,jj) = xx + dt1*scale_rot* (+lw(1)*yy+lw(2)*zz)
552 eigenvec(2,jj) = yy + dt1*scale_rot* (-lw(1)*xx+lw(3)*zz)
553 eigenvec(3,jj) = zz + dt1*scale_rot* (-lw(2)*xx-lw(3)*yy)
554 ENDDO
555 ! save it in data structure (RESTART)
556 ale%GRID%flow_tracking_data%EIGENVEC(1:3,1)=eigenvec(1:3,1)
557 ale%GRID%flow_tracking_data%EIGENVEC(1:3,2)=eigenvec(1:3,2)
558 ale%GRID%flow_tracking_data%EIGENVEC(1:3,3)=eigenvec(1:3,3)
559 ENDIF
560
561C-----------------------------------------------
562 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215
ale_mod
Definition ale_mod.F:149
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:249
my_barrier
subroutine my_barrier
Definition machine.F:31