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