OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_muscl_gradients.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| multi_muscl_gradients ../engine/source/multifluid/multi_muscl_gradients.F
25!||--- called by ------------------------------------------------------
26!|| multi_timeevolution ../engine/source/multifluid/multi_timeevolution.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| limiter ../engine/source/multifluid/multi_muscl_gradients.F
30!|| startime ../engine/source/system/timer_mod.F90
31!|| stoptime ../engine/source/system/timer_mod.F90
32!||--- uses -----------------------------------------------------
33!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!|| element_mod ../common_source/modules/elements/element_mod.F90
36!|| initbuf_mod ../engine/share/resol/initbuf.F
37!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
38!|| timer_mod ../engine/source/system/timer_mod.F90
39!||====================================================================
40 SUBROUTINE multi_muscl_gradients(TIMERS,ELBUF_TAB, IPARG, ITASK, IXS, IXQ, IXTG,
41 . PM, IPM, MULTI_FVM, ALE_CONNECTIVITY, WGRID, XGRID, ITAB, NBMAT, CURRENT_TIME, BUFMAT)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE timer_mod
46 USE initbuf_mod
47 USE elbufdef_mod
48 USE multi_fvm_mod
49 USE ale_connectivity_mod
50 use element_mod , only : nixs,nixq,nixtg
51C-----------------------------------------------
52C I m p l i c i t T y p e s
53C-----------------------------------------------
54#include "implicit_f.inc"
55C-----------------------------------------------
56C C o m m o n B l o c k s
57C-----------------------------------------------
58#include "com01_c.inc"
59! NUMELS
60#include "param_c.inc"
61#include "task_c.inc"
62! DTFAC1
63#include "macro.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 TYPE(timer_), INTENT(INOUT) :: TIMERS
68 TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
69 INTEGER, INTENT(IN) :: IPARG(NPARG, *)
70 INTEGER, INTENT(IN) :: ITASK ! SMP TASK
71 INTEGER, INTENT(IN) :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
72 INTEGER, INTENT(IN) :: IPM(NPROPMI, *)
73 my_real, INTENT(IN) :: pm(npropm, *)
74 TYPE(multi_fvm_struct), INTENT(INOUT) :: MULTI_FVM
75 my_real, INTENT(IN) :: wgrid(3, *), xgrid(3, *), current_time
76 INTEGER, INTENT(IN) :: ITAB(*), NBMAT
77 my_real, INTENT(INOUT) :: bufmat(*)
78 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECTIVITY
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER :: II, I, J, NB_FACE, NB_FACE2, KFACE, IFACE, NG, IMAT
83 my_real :: valk, vall, nx(6), ny(6), nz(6),
84 . xk(3), xf(3, 6), delta, lim_rho, xl(3, 6)
85 my_real :: maxval_rho, minval_rho,
86 . maxval_u, minval_u,
87 . maxval_v, minval_v,
88 . maxval_w, minval_w, lim_velx, lim_vely, lim_velz, lim_vel
89 my_real :: maxval_pres, minval_pres, lim_pres
90 my_real :: phase_maxval_rho(nbmat), phase_minval_rho(nbmat),
91 . phase_maxval_alpha(nbmat), phase_minval_alpha(nbmat),
92 . phase_maxval_pres(nbmat), phase_minval_pres(nbmat),
93 . phase_lim_alpha(nbmat), phase_lim_rho(nbmat), phase_lim_pres(nbmat)
94 my_real :: peps, meps, limface(6)
95 TYPE(g_bufel_), POINTER :: GBUF
96 TYPE(L_BUFEL_), POINTER :: LBUF
97 INTEGER :: VOIS(6), FACE_LIST(6), NEL, ITY, ISOLNOD, NFT, MATLAW,
98 . node1, node2, node3, node4
99 my_real :: sumx, sumy, sumz, vol, one_over_vol, lim_alpha
100 my_real :: eps_int
101 my_real, DIMENSION(:), POINTER :: volume
102 my_real :: x1(3), x2(3), x3(3), x4(3)
103 my_real :: mat(3, 3), rhs(3), sol(3), matm1(3, 3), det, mattest(3, 3)
104 INTEGER :: IAD, LGTH
105
106
107
108 IF(itask == 0) CALL startime(timers,macro_timer_muscl)
109 eps_int = five * em02
110
111 DO ng = itask + 1, ngroup, nthread
112 matlaw = iparg(1, ng)
113 IF (matlaw == 151) THEN
114 nel = iparg(2, ng)
115 nft = iparg(3, ng)
116 ity = iparg(5, ng)
117 isolnod = iparg(28, ng)
118
119 gbuf => elbuf_tab(ng)%GBUF
120 IF (multi_fvm%SYM == 0) THEN
121 volume => gbuf%VOL(1:nel)
122 ELSE
123 volume => gbuf%AREA(1:nel)
124 ENDIF
125
126 peps = em10
127 meps = -em10
128
129 DO iface = 1, 6
130 face_list(iface) = iface
131 ENDDO
132
133 nb_face = -1
134 nb_face2 = -1
135 SELECT CASE (multi_fvm%SYM)
136 CASE (0)
137C Solids
138 nb_face = 6
139 nb_face2 = 6
140 IF (isolnod == 4) THEN
141C Tetra
142 nb_face = 4
143 face_list(1) = 5
144 face_list(2) = 6
145 face_list(3) = 2
146 face_list(4) = 4
147 face_list(5) = -1
148 face_list(6) = -1
149 ENDIF
150 CASE (1, 2)
151 IF (ity == 2) THEN
152C QUADS
153 nb_face = 4
154 nb_face2 = 4
155 ELSEIF (ity == 7) THEN
156C TRIANGLES
157 nb_face = 3
158 nb_face2 = 3
159 ENDIF
160 CASE DEFAULT
161 CALL arret(2)
162 END SELECT
163
164 DO ii = 1, nel
165 i = ii + nft
166 iad = ale_connectivity%ee_connect%iad_connect(i)
167! NB_FACE = ALE_CONNECTIVITY%ee_connect%iad_connect(I+1) -
168! . ALE_CONNECTIVITY%ee_connect%iad_connect(I)
169 vois(1:6) = -1
170 nx(1:6) = -ep30
171 ny(1:6) = -ep30
172 nz(1:6) = -ep30
173 vol = volume(ii)
174 one_over_vol = one / vol
175C centroid element
176 xk(1:3) = multi_fvm%ELEM_DATA%CENTROID(1:3, i)
177C Matrice
178 mat(1:3, 1:3) = zero
179 DO iface = 1, nb_face
180 kface = face_list(iface)
181 j = ale_connectivity%ee_connect%connected(iad + kface - 1)
182 xf(1:3, kface) = multi_fvm%FACE_DATA%CENTROID(1:3, kface, i)
183 IF (j > 0) THEN
184 vois(kface) = j
185 xl(1:3, kface) = multi_fvm%ELEM_DATA%CENTROID(1:3, j)
186 ELSE
187 xl(1:3, kface) = two * xf(1:3, kface) - xk(1:3)
188 ENDIF
189 nx(kface) = multi_fvm%FACE_DATA%NORMAL(1, kface, i)
190 ny(kface) = multi_fvm%FACE_DATA%NORMAL(2, kface, i)
191 nz(kface) = multi_fvm%FACE_DATA%NORMAL(3, kface, i)
192 mat(1, 1) = mat(1, 1) + (xl(1, kface) - xk(1)) * (xl(1, kface) - xk(1))
193 mat(1, 2) = mat(1, 2) + (xl(1, kface) - xk(1)) * (xl(2, kface) - xk(2))
194 mat(1, 3) = mat(1, 3) + (xl(1, kface) - xk(1)) * (xl(3, kface) - xk(3))
195 mat(2, 1) = mat(2, 1) + (xl(2, kface) - xk(2)) * (xl(1, kface) - xk(1))
196 mat(2, 2) = mat(2, 2) + (xl(2, kface) - xk(2)) * (xl(2, kface) - xk(2))
197 mat(2, 3) = mat(2, 3) + (xl(2, kface) - xk(2)) * (xl(3, kface) - xk(3))
198 mat(3, 1) = mat(3, 1) + (xl(3, kface) - xk(3)) * (xl(1, kface) - xk(1))
199 mat(3, 2) = mat(3, 2) + (xl(3, kface) - xk(3)) * (xl(2, kface) - xk(2))
200 mat(3, 3) = mat(3, 3) + (xl(3, kface) - xk(3)) * (xl(3, kface) - xk(3))
201 ENDDO
202
203 IF (multi_fvm%SYM /= 0) THEN
204 mat(1, 1) = one
205 mat(1, 2) = zero
206 mat(1, 3) = zero
207 mat(2, 1) = zero
208 mat(3, 1) = zero
209 ENDIF
210
211 det = mat(1, 1) * mat(2, 2) * mat(3, 3) - mat(1, 1) * mat(2, 3)**2
212 . - mat(2, 2) * mat(1, 3)**2 - mat(3, 3) * mat(1, 2)**2
213 . + two * mat(1, 2) * mat(1, 3) * mat(2, 3)
214 IF (det /= zero) THEN
215 det = one / det
216 ELSE
217 CALL arret(2)
218 ENDIF
219
220 matm1(1, 1) = mat(2, 2) * mat(3, 3) - mat(2, 3)**2
221 matm1(1, 2) = -mat(1, 2) * mat(3, 3) + mat(1, 3) * mat(2, 3)
222 matm1(1, 3) = mat(1, 2) * mat(2, 3) - mat(1, 3) * mat(2, 2)
223 matm1(2, 2) = mat(1, 1) * mat(3, 3) - mat(1, 3)**2
224 matm1(2, 3) = -mat(1, 1) * mat(2, 3) + mat(1, 2) * mat(1, 3)
225 matm1(3, 3) = mat(1, 1) * mat(2, 2) - mat(1, 2)**2
226 matm1(2, 1) = matm1(1, 2)
227 matm1(3, 1) = matm1(1, 3)
228 matm1(3, 2) = matm1(2, 3)
229
230 matm1(1:3, 1:3) = matm1(1:3, 1:3) * det
231
232C MATTEST(1:3, 1:3) = MATMUL(MAT(1:3, 1:3), MATM1(1:3, 1:3))
233
234 IF (multi_fvm%SYM == 1) THEN
235C Use planar surface here
236 IF (ity == 2) THEN
237C QUADS
238 node1 = ixq(2, i)
239 node2 = ixq(3, i)
240 node3 = ixq(4, i)
241 node4 = ixq(5, i)
242 x1(1:3) = xgrid(1:3, node1)
243 x2(1:3) = xgrid(1:3, node2)
244 x3(1:3) = xgrid(1:3, node3)
245 x4(1:3) = xgrid(1:3, node4)
246 ELSE IF (ity == 7) THEN
247C TRIANGLES
248 node1 = ixtg(2, i)
249 node2 = ixtg(3, i)
250 node3 = ixtg(4, i)
251 x1(1:3) = xgrid(1:3, node1)
252 x2(1:3) = xgrid(1:3, node2)
253 x3(1:3) = xgrid(1:3, node3)
254 ENDIF
255 ENDIF
256
257C Velocity gradient
258 maxval_u = multi_fvm%VEL(1, i)
259 minval_u = maxval_u
260 maxval_v = multi_fvm%VEL(2, i)
261 minval_v = maxval_v
262 maxval_w = multi_fvm%VEL(3, i)
263 minval_w = maxval_w
264
265 IF (nbmat == 1) THEN
266C =========
267C MONOFLUID
268C =========
269C Density gradient
270 maxval_rho = multi_fvm%RHO(i)
271 minval_rho = maxval_rho
272C Gradient press
273 maxval_pres = multi_fvm%EINT(i)
274 minval_pres = maxval_pres
275 ELSE
276C ==========
277C MULTIFLUID
278C ==========
279 DO imat = 1, nbmat
280C Density gradient
281 phase_maxval_rho(imat) = multi_fvm%PHASE_RHO(imat, i)
282 phase_minval_rho(imat) = phase_maxval_rho(imat)
283C Gradient press
284 phase_maxval_pres(imat) = multi_fvm%PHASE_EINT(imat, i)
285 phase_minval_pres(imat) = phase_maxval_pres(imat)
286C Volume fraction gradient
287 phase_maxval_alpha(imat) = multi_fvm%PHASE_ALPHA(imat, i)
288 phase_minval_alpha(imat) = phase_maxval_alpha(imat)
289 ENDDO
290 ENDIF
291
292 IF (multi_fvm%MUSCL == 1) THEN
293C Velocity gradient -> X component
294 valk = multi_fvm%VEL(1, i)
295 rhs(1:3) = zero
296 DO iface = 1, nb_face
297 kface = face_list(iface)
298 j = vois(kface)
299 vall = valk
300 IF (j > 0) THEN
301 vall = multi_fvm%VEL(1, j)
302 maxval_u = max(maxval_u, vall)
303 minval_u = min(minval_u, vall)
304 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
305 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
306 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
307 ENDIF
308 ENDDO
309
310 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
311 sumx = -sol(1)
312 sumy = -sol(2)
313 sumz = -sol(3)
314C Component X limitation
315 limface(1:6) = one
316 DO iface = 1, nb_face
317 kface = face_list(iface)
318 delta = sumx * (xf(1, kface) - xk(1)) +
319 . sumy * (xf(2, kface) - xk(2)) +
320 . sumz * (xf(3, kface) - xk(3))
321 IF (delta > zero) THEN
322 limface(kface) = half * (maxval_u - valk) / delta
323 ELSE IF (delta < zero) THEN
324 limface(kface) = half * (minval_u - valk) / delta
325 ENDIF
326 CALL limiter(limface(kface), one)
327 ENDDO
328 lim_velx = minval(limface(1:6))
329 multi_fvm%GRAD_U(1, i) = sumx * lim_velx
330 multi_fvm%GRAD_U(2, i) = sumy * lim_velx
331 multi_fvm%GRAD_U(3, i) = sumz * lim_velx
332C Velocity gradient -> Y component
333 valk = multi_fvm%VEL(2, i)
334 rhs(1:3) = zero
335 DO iface = 1, nb_face
336 kface = face_list(iface)
337 j = vois(kface)
338 vall = valk
339 IF (j > 0) THEN
340 vall = multi_fvm%VEL(2, j)
341 maxval_v = max(maxval_v, vall)
342 minval_v = min(minval_v, vall)
343 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
344 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
345 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
346 ENDIF
347 ENDDO
348 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
349 sumx = -sol(1)
350 sumy = -sol(2)
351 sumz = -sol(3)
352C Component Y limitation
353 limface(1:6) = one
354 DO iface = 1, nb_face
355 kface = face_list(iface)
356 delta = sumx * (xf(1, kface) - xk(1)) +
357 . sumy * (xf(2, kface) - xk(2)) +
358 . sumz * (xf(3, kface) - xk(3))
359 IF (delta > zero) THEN
360 limface(kface) = half * (maxval_v - valk) / delta
361 ELSE IF (delta < zero) THEN
362 limface(kface) = half * (minval_v - valk) / delta
363 ENDIF
364 CALL limiter(limface(kface), one)
365 ENDDO
366 lim_vely = minval(limface(1:6))
367 multi_fvm%GRAD_V(1, i) = sumx * lim_vely
368 multi_fvm%GRAD_V(2, i) = sumy * lim_vely
369 multi_fvm%GRAD_V(3, i) = sumz * lim_vely
370C Velocity gradient -> Z component
371 valk = multi_fvm%VEL(3, i)
372 rhs(1:3) = zero
373 DO iface = 1, nb_face
374 kface = face_list(iface)
375 j = vois(kface)
376 vall = valk
377 IF (j > 0) THEN
378 vall = multi_fvm%VEL(3, j)
379 maxval_w = max(maxval_w, vall)
380 minval_w = min(minval_w, vall)
381 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
382 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
383 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
384 ENDIF
385 ENDDO
386 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
387 sumx = -sol(1)
388 sumy = -sol(2)
389 sumz = -sol(3)
390C Component Z limitation
391 limface(1:6) = one
392 DO iface = 1, nb_face
393 kface = face_list(iface)
394 delta = sumx * (xf(1, kface) - xk(1)) +
395 . sumy * (xf(2, kface) - xk(2)) +
396 . sumz * (xf(3, kface) - xk(3))
397 IF (delta > zero) THEN
398 limface(kface) = half * (maxval_w - valk) / delta
399 ELSE IF (delta < zero) THEN
400 limface(kface) = half * (minval_w - valk) / delta
401 ENDIF
402 CALL limiter(limface(kface), one)
403 ENDDO
404 lim_velz = minval(limface(1:6))
405 multi_fvm%GRAD_W(1, i) = sumx * lim_velz
406 multi_fvm%GRAD_W(2, i) = sumy * lim_velz
407 multi_fvm%GRAD_W(3, i) = sumz * lim_velz
408 ENDIF !MUSCL == 1
409
410 IF (nbmat == 1) THEN
411C ==========
412C MONOFLUID
413C ==========
414C Density gradient
415 IF (multi_fvm%MUSCL == 1) THEN
416 valk = multi_fvm%RHO(i)
417 rhs(1:3) = zero
418 DO iface = 1, nb_face
419 kface = face_list(iface)
420 j = vois(kface)
421 vall = valk
422 IF (j > 0) THEN
423 vall = multi_fvm%RHO(j)
424 maxval_rho = max(maxval_rho, vall)
425 minval_rho = min(minval_rho, vall)
426 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
427 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
428 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
429 ENDIF
430 ENDDO
431 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
432 sumx = -sol(1)
433 sumy = -sol(2)
434 sumz = -sol(3)
435C Density gradient limitation
436 limface(1:6) = one
437 DO iface = 1, nb_face
438 kface = face_list(iface)
439 delta = sumx * (xf(1, kface) - xk(1)) +
440 . sumy * (xf(2, kface) - xk(2)) +
441 . sumz * (xf(3, kface) - xk(3))
442 IF (delta > zero) THEN
443 limface(kface) = half * (maxval_rho - valk) / delta
444 ELSE IF (delta < zero) THEN
445 limface(kface) = half * (minval_rho - valk) / delta
446 ENDIF
447 CALL limiter(limface(kface), one)
448 ENDDO
449 lim_rho = minval(limface(1:6))
450 multi_fvm%GRAD_RHO(1, i) = sumx * lim_rho
451 multi_fvm%GRAD_RHO(2, i) = sumy * lim_rho
452 multi_fvm%GRAD_RHO(3, i) = sumz * lim_rho
453C Gradient press
454 valk = multi_fvm%EINT(i)
455 rhs(1:3) = zero
456 DO iface = 1, nb_face
457 kface = face_list(iface)
458 j = vois(kface)
459 vall = valk
460 IF (j > 0) THEN
461 vall = multi_fvm%EINT(j)
462 maxval_pres = max(maxval_pres, vall)
463 minval_pres = min(minval_pres, vall)
464 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
465 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
466 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
467 ENDIF
468 ENDDO
469 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
470 sumx = -sol(1)
471 sumy = -sol(2)
472 sumz = -sol(3)
473
474C pressure limitation
475 limface(1:6) = one
476 DO iface = 1, nb_face
477 kface = face_list(iface)
478 delta = sumx * (xf(1, kface) - xk(1)) +
479 . sumy * (xf(2, kface) - xk(2)) +
480 . sumz * (xf(3, kface) - xk(3))
481 IF (delta > zero) THEN
482 limface(kface) = half * (maxval_pres - valk) / delta
483 ELSE IF (delta < zero) THEN
484 limface(kface) = half * (minval_pres - valk) / delta
485 ENDIF
486 CALL limiter(limface(kface), one)
487 ENDDO
488 lim_pres = minval(limface(1:6))
489 multi_fvm%GRAD_PRES(1, i) = sumx * lim_pres
490 multi_fvm%GRAD_PRES(2, i) = sumy * lim_pres
491 multi_fvm%GRAD_PRES(3, i) = sumz * lim_pres
492 ENDIF ! MUSCL == 1
493 ELSE
494C ==========
495C MULTIFLUID
496C ==========
497 DO imat = 1, nbmat
498C Volume fraction gradient
499 valk = multi_fvm%PHASE_ALPHA(imat, i)
500 rhs(1:3) = zero
501 DO iface = 1, nb_face
502 kface = face_list(iface)
503 j = vois(kface)
504 vall = valk
505 IF (j > 0) THEN
506 vall = multi_fvm%PHASE_ALPHA(imat, j)
507 phase_maxval_alpha(imat) = max(phase_maxval_alpha(imat), vall)
508 phase_minval_alpha(imat) = min(phase_minval_alpha(imat), vall)
509 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
510 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
511 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
512 ENDIF
513 ENDDO
514 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
515 sumx = -sol(1)
516 sumy = -sol(2)
517 sumz = -sol(3)
518C Volume fraction gradient limitation
519 limface(1:6) = one
520 DO iface = 1, nb_face
521 kface = face_list(iface)
522 delta = sumx * (xf(1, kface) - xk(1)) +
523 . sumy * (xf(2, kface) - xk(2)) +
524 . sumz * (xf(3, kface) - xk(3))
525 IF (delta > zero) THEN
526 limface(kface) = half * (phase_maxval_alpha(imat) - valk) / delta
527 ELSE IF (delta < zero) THEN
528 limface(kface) = half * (phase_minval_alpha(imat) - valk) / delta
529 ENDIF
530 CALL limiter(limface(kface), multi_fvm%BETA)
531 ENDDO
532 phase_lim_alpha(imat) = minval(limface(1:6))
533 multi_fvm%PHASE_GRAD_ALPHA(1, imat, i) = sumx * phase_lim_alpha(imat)
534 multi_fvm%PHASE_GRAD_ALPHA(2, imat, i) = sumy * phase_lim_alpha(imat)
535 multi_fvm%PHASE_GRAD_ALPHA(3, imat, i) = sumz * phase_lim_alpha(imat)
536C Density gradient
537 IF (multi_fvm%MUSCL == 1) THEN
538 valk = multi_fvm%PHASE_RHO(imat, i)
539 rhs(1:3) = zero
540 DO iface = 1, nb_face
541 kface = face_list(iface)
542 j = vois(kface)
543 vall = valk
544 IF (j > 0) THEN
545 vall = multi_fvm%PHASE_RHO(imat, j)
546 phase_maxval_rho(imat) = max(phase_maxval_rho(imat), vall)
547 phase_minval_rho(imat) = min(phase_minval_rho(imat), vall)
548 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
549 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
550 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
551 ENDIF
552 ENDDO
553 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
554 sumx = -sol(1)
555 sumy = -sol(2)
556 sumz = -sol(3)
557C Density gradient limitation
558 limface(1:6) = one
559 DO iface = 1, nb_face
560 kface = face_list(iface)
561 delta = sumx * (xf(1, kface) - xk(1)) +
562 . sumy * (xf(2, kface) - xk(2)) +
563 . sumz * (xf(3, kface) - xk(3))
564 IF (delta > zero) THEN
565 limface(kface) = half * (phase_maxval_rho(imat) - valk) / delta
566 ELSE IF (delta < zero) THEN
567 limface(kface) = half * (phase_minval_rho(imat) - valk) / delta
568 ENDIF
569 IF (valk + delta <= zero) THEN
570 limface(kface) = zero
571 ELSE
572 CALL limiter(limface(kface), one)
573 ENDIF
574 ENDDO
575 phase_lim_rho(imat) = minval(limface(1:6))
576 multi_fvm%PHASE_GRAD_RHO(1, imat, i) = sumx * phase_lim_rho(imat)
577 multi_fvm%PHASE_GRAD_RHO(2, imat, i) = sumy * phase_lim_rho(imat)
578 multi_fvm%PHASE_GRAD_RHO(3, imat, i) = sumz * phase_lim_rho(imat)
579C Gradient press
580 valk = multi_fvm%PHASE_EINT(imat, i)
581 rhs(1:3) = zero
582 DO iface = 1, nb_face
583 kface = face_list(iface)
584 j = vois(kface)
585 vall = valk
586 IF (j > 0) THEN
587 vall = multi_fvm%PHASE_EINT(imat, j)
588 phase_maxval_pres(imat) = max(phase_maxval_pres(imat), vall)
589 phase_minval_pres(imat) = min(phase_minval_pres(imat), vall)
590 rhs(1) = rhs(1) + (valk - vall) * (xl(1, kface) - xk(1))
591 rhs(2) = rhs(2) + (valk - vall) * (xl(2, kface) - xk(2))
592 rhs(3) = rhs(3) + (valk - vall) * (xl(3, kface) - xk(3))
593 ENDIF
594 ENDDO
595 sol(1:3) = matmul(matm1(1:3, 1:3), rhs(1:3))
596 sumx = -sol(1)
597 sumy = -sol(2)
598 sumz = -sol(3)
599C Gradient Limitation Pressure
600 limface(1:6) = one
601 DO iface = 1, nb_face
602 kface = face_list(iface)
603 delta = sumx * (xf(1, kface) - xk(1)) +
604 . sumy * (xf(2, kface) - xk(2)) +
605 . sumz * (xf(3, kface) - xk(3))
606 IF (delta > zero) THEN
607 limface(kface) = half * (phase_maxval_pres(imat) - valk) / delta
608 ELSE IF (delta < zero) THEN
609 limface(kface) = half * (phase_minval_pres(imat) - valk) / delta
610 ENDIF
611 CALL limiter(limface(kface), one)
612 ENDDO
613 phase_lim_pres(imat) = minval(limface(1:6))
614 multi_fvm%PHASE_GRAD_PRES(1, imat, i) = sumx * phase_lim_pres(imat)
615 multi_fvm%PHASE_GRAD_PRES(2, imat, i) = sumy * phase_lim_pres(imat)
616 multi_fvm%PHASE_GRAD_PRES(3, imat, i) = sumz * phase_lim_pres(imat)
617 ENDIF ! MUSCL == 1
618 ENDDO ! IMAT
619 ENDIF ! NBMAT
620 ENDDO !II = 1, NEL
621 ENDIF !MATLAW == 151
622 ENDDO !NG
623 IF(itask == 0) CALL stoptime(timers,macro_timer_muscl)
624
625 END SUBROUTINE multi_muscl_gradients
626
627!||====================================================================
628!|| limiter ../engine/source/multifluid/multi_muscl_gradients.F
629!||--- called by ------------------------------------------------------
630!|| multi_muscl_gradients ../engine/source/multifluid/multi_muscl_gradients.F
631!||====================================================================
632 SUBROUTINE limiter(PHI, BETA)
633#include "implicit_f.inc"
634 my_real, INTENT(INOUT) :: phi
635 my_real, INTENT(IN) :: beta
636 my_real :: lim
637
638 lim = max(zero, max(min(one, beta * phi), min(phi, beta)))
639 phi = lim
640 END SUBROUTINE limiter
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
multi_muscl_gradients
subroutine multi_muscl_gradients(timers, elbuf_tab, iparg, itask, ixs, ixq, ixtg, pm, ipm, multi_fvm, ale_connectivity, wgrid, xgrid, itab, nbmat, current_time, bufmat)
Definition multi_muscl_gradients.F:42
limiter
subroutine limiter(phi, beta)
Definition multi_muscl_gradients.F:633
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
initbuf_mod
Definition initbuf.F:160
arret
subroutine arret(nn)
Definition arret.F:86
startime
subroutine startime(event, itask)
Definition timer.F:93
stoptime
subroutine stoptime(event, itask)
Definition timer.F:135
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256