OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_pressure_equilibrium.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_pressure_equilibrium ../engine/source/multifluid/multi_pressure_equilibrium.F
25!||--- called by ------------------------------------------------------
26!|| alemain ../engine/source/ale/alemain.F
27!||--- calls -----------------------------------------------------
28!|| multi_submatlaw ../engine/source/multifluid/multi_submatlaw.F
29!||--- uses -----------------------------------------------------
30!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
31!|| initbuf_mod ../engine/share/resol/initbuf.F
32!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
33!|| multi_fvm_mod ../common_source/modules/ale/multi_fvm_mod.F90
34!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.F
35!||====================================================================
36 SUBROUTINE multi_pressure_equilibrium(TIMESTEP, ELBUF_TAB, IPARG, ITASK, IXS, IXQ, IXTG,
37 . PM, IPM, MULTI_FVM, CURRENT_TIME, BUFMAT, NPF, TF, NUMMAT,MATPARAM)
38C-----------------------------------------------
39C M o d u l e s
40C-----------------------------------------------
41 USE initbuf_mod
42 USE elbufdef_mod
43 USE multi_fvm_mod
44 USE matparam_def_mod, ONLY : matparam_struct_
45 USE multi_submatlaw_mod , ONLY : multi_submatlaw
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "com01_c.inc"
54#include "param_c.inc"
55#include "task_c.inc"
56#include "mvsiz_p.inc"
57#include "vect01_c.inc"
58#include "tabsiz_c.inc"
59C-----------------------------------------------
60C D u m m y A r g u m e n t s
61C-----------------------------------------------
62 TYPE(matparam_struct_),DIMENSION(NUMMAT),INTENT(IN) :: MATPARAM !< material buffer
63 INTEGER,INTENT(IN) :: NUMMAT !< number of material law (array size for IPM and PM)
64 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
65 INTEGER, INTENT(IN) :: IPARG(NPARG, NGROUP)
66 INTEGER, INTENT(IN) :: ITASK ! SMP TASK
67 INTEGER, INTENT(IN), TARGET :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
68 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT)
69 my_real, INTENT(IN) :: pm(npropm, nummat)
70 TYPE(multi_fvm_struct), INTENT(INOUT), TARGET :: MULTI_FVM
71 my_real, INTENT(IN) :: current_time, timestep
72 my_real, INTENT(INOUT) :: bufmat(*)
73 INTEGER,INTENT(IN)::NPF(SNPC)
74 my_real,INTENT(IN)::tf(stf)
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 TYPE(g_bufel_), POINTER :: GBUF
79 TYPE(l_bufel_), POINTER :: LBUF
80 TYPE(buf_eos_), POINTER :: EBUF
81 INTEGER, DIMENSION(:, :), POINTER :: IX
82 INTEGER :: NG, II, I
83 INTEGER :: NEL
84 INTEGER :: NBMAT, IMAT, LOCAL_MATID, MATLAW
85 my_real :: vol_frac
86 my_real, DIMENSION(:), POINTER :: rho, ssp, vol, pres, eint, tau, speeint, feint, fp
87 INTEGER :: ITER, MAX_ITER, REMAINING_ELTS
88 my_real :: tol, tole
89 LOGICAL :: CONT
90 my_real :: grun(mvsiz, multi_fvm%NBMAT), coef(mvsiz)
91 my_real :: fe(mvsiz)
92 INTEGER :: NBMAT_INCELL(MVSIZ)
93 my_real :: one_over_grun
94 my_real :: error(mvsiz), relaxed_pressure(mvsiz)
95 my_real :: pres_incr(mvsiz), tau_incr(mvsiz, multi_fvm%NBMAT)
96 INTEGER :: CV_INDICATOR(MVSIZ)
97 my_real, TARGET :: volwt(mvsiz, multi_fvm%NBMAT), eintwt(mvsiz, multi_fvm%NBMAT),
98 . rhowt(mvsiz, multi_fvm%NBMAT), preswt(mvsiz, multi_fvm%NBMAT),
99 . sspwt(mvsiz, multi_fvm%NBMAT),
100 . masswt(mvsiz, multi_fvm%NBMAT), tauwt(mvsiz, multi_fvm%NBMAT),
101 . speeintwt(mvsiz, multi_fvm%NBMAT), alphawt(mvsiz, multi_fvm%NBMAT),
102 . mfracwt(mvsiz, multi_fvm%NBMAT),
103 . fpwt(mvsiz, multi_fvm%NBMAT), fewt(mvsiz, multi_fvm%NBMAT),
104 . fe0wt(mvsiz, multi_fvm%NBMAT), tau0wt(mvsiz, multi_fvm%NBMAT)
105 my_real :: dpdtau, dpde, p0(mvsiz), coef1, lim_tau(mvsiz, multi_fvm%NBMAT), lim
106 INTEGER :: NIX,NVAREOS
107 LOGICAL :: IDBG
108 INTEGER :: NVARTMP_EOS
109 my_real :: pshift
110 my_real :: sigold(mvsiz*6)
111C-----------------------------------------------
112C B e g i n n i n g o f s u b r o u t i n e
113C-----------------------------------------------
114 idbg = .false.
115 NULLIFY(ix)
116
117 pshift = multi_fvm%PRES_SHIFT
118
119 DO ng = itask + 1, ngroup, nthread
120 mtn = iparg(1, ng)
121 IF (mtn == 151) THEN
122 nel = iparg(2, ng)
123 nft = iparg(3, ng)
124 ity = iparg(5, ng)
125 lft = 1
126 llt = nel
127 IF (multi_fvm%SYM == 0) THEN
128 ix => ixs(1:nixs, 1 + nft:nel + nft)
129 nix = nixs
130 ELSEIF (ity == 2) THEN
131C QUADS
132 ix => ixq(1:nixq, 1 + nft:nel + nft)
133 nix = nixq
134 ELSE !IF (ITY == 7) THEN
135C TRIANGLES
136 ix => ixtg(1:nixtg, 1 + nft:nel + nft)
137 nix = nixtg
138 ENDIF
139C Multifluid law number, get the number of materials
140 nbmat = multi_fvm%NBMAT
141C Global buffer of the current group
142 gbuf => elbuf_tab(ng)%GBUF
143 sigold(1 + 0 * nel : 1 * nel)=gbuf%SIG(1 + 0 * nel : 1 * nel)
144 sigold(1 + 1 * nel : 2 * nel)=gbuf%SIG(1 + 1 * nel : 2 * nel)
145 sigold(1 + 2 * nel : 3 * nel)=gbuf%SIG(1 + 2 * nel : 3 * nel)
146 sigold(1 + 3 * nel : 4 * nel)=gbuf%SIG(1 + 3 * nel : 4 * nel)
147 sigold(1 + 4 * nel : 5 * nel)=gbuf%SIG(1 + 4 * nel : 5 * nel)
148 sigold(1 + 5 * nel : 6 * nel)=gbuf%SIG(1 + 5 * nel : 6 * nel)
149 IF (nbmat == 1) THEN
150 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1, 1, 1)
151 imat = 1
152 IF(ASSOCIATED(ix)) THEN
153 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
154 ELSE
155 local_matid = 1
156 ENDIF
157 matlaw = matparam(local_matid)%ILAW
158 pres => multi_fvm%PRES(1 + nft : nel + nft)
159 ssp => multi_fvm%SOUND_SPEED(1 + nft : nel + nft)
160 vol => gbuf%VOL
161 eint => multi_fvm%EINT(1 + nft : nel + nft)
162 rho => multi_fvm%RHO(1 + nft : nel + nft)
163 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1, 1, 1)
164 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
165 nvartmp_eos = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS
166 CALL multi_submatlaw(
167 1 1, matlaw, local_matid, nel,
168 2 eint, pres, rho, ssp,
169 3 vol, grun(1:nel,1),pm, ipm,
170 4 npropm, npropmi, bufmat, gbuf%OFF,
171 5 gbuf%TEMP, lbuf%BFRAC, gbuf%TB, gbuf%DELTAX,
172 6 current_time, sigold, snpc, stf,
173 7 npf, tf, ebuf%VAR, nvareos,
174 7 matparam(local_matid), nvartmp_eos, ebuf%VARTMP, nummat ,
175 8 gbuf%ABURN)
176 DO ii = 1, nel
177 IF (multi_fvm%SOUND_SPEED(ii + nft) > zero) THEN
178 multi_fvm%SOUND_SPEED(ii + nft) = sqrt(multi_fvm%SOUND_SPEED(ii + nft))
179 ELSE
180 multi_fvm%SOUND_SPEED(ii + nft) = em10
181 ENDIF
182 ENDDO
183 IF (matlaw == 5) THEN
184 DO ii = 1, nel
185 i = ii + nft
186 multi_fvm%BFRAC(1, i) = lbuf%BFRAC(ii)
187 ENDDO
188 ENDIF
189 gbuf%SIG(1 + 0 * nel : 1 * nel) = -pres(1:nel)
190 gbuf%SIG(1 + 1 * nel : 2 * nel) = gbuf%SIG(1 + 0 * nel : 1 * nel)
191 gbuf%SIG(1 + 2 * nel : 3 * nel) = gbuf%SIG(1 + 0 * nel : 1 * nel)
192 ELSE
193C=======================================================================
194C Count number of mat present in cell
195C========================
196 nbmat_incell(1:nel) = 0
197 DO imat = 1, nbmat
198 DO ii = 1, nel
199 i = ii + nft
200 IF (multi_fvm%PHASE_ALPHA(imat, i) > zero .AND.
201 . multi_fvm%PHASE_ALPHA(imat, i) * multi_fvm%PHASE_RHO(imat, i)
202 . / multi_fvm%RHO(i) > em14) THEN
203 ! At this stage, only significant materials have a
204 ! positive volume. Filtering has been made in MULTI_EVOLVE_PARTIAL
205 nbmat_incell(ii) = nbmat_incell(ii) + 1
206 ENDIF
207 ENDDO
208 ENDDO
209C For cells containing only one material, set eint to global eint
210 DO imat = 1, nbmat
211 DO ii = 1, nel
212 i = ii + nft
213 IF (nbmat_incell(ii) == 1 .AND.
214 . multi_fvm%PHASE_ALPHA(imat, i) > zero) THEN
215 multi_fvm%PHASE_EINT(imat, i) = multi_fvm%EINT(i)
216 ENDIF
217 ENDDO
218 ENDDO
219C=======================================================================
220C PRESSURE RELAXATION
221C=======================================================================
222C As temporary variables, in order not to lose values in case of convergence failure
223 DO imat = 1, nbmat
224 DO ii = 1, nel
225 i = ii + nft
226C Volume
227 volwt(ii, imat) = gbuf%VOL(ii) * multi_fvm%PHASE_ALPHA(imat, i)
228C Volume fraction
229 alphawt(ii, imat) = multi_fvm%PHASE_ALPHA(imat, i)
230C Internal energy
231 eintwt(ii, imat) = multi_fvm%PHASE_EINT(imat, i)
232
233
234C Density
235 rhowt(ii, imat) = multi_fvm%PHASE_RHO(imat, i)
236C Mass fraction
237 mfracwt(ii, imat) = alphawt(ii, imat) * rhowt(ii, imat) / gbuf%RHO(ii)
238C Specific volume
239 tauwt(ii, imat) = zero
240C Specific internal energy
241 speeintwt(ii, imat) = zero
242 IF (rhowt(ii, imat) > zero) THEN
243 tauwt(ii, imat) = one / rhowt(ii, imat)
244 speeintwt(ii, imat) = eintwt(ii, imat) / rhowt(ii, imat)
245 ENDIF
246 tau0wt(ii, imat) = tauwt(ii, imat)
247C Normalized mass (alpha rho)
248 masswt(ii, imat) = rhowt(ii, imat) * alphawt(ii, imat)
249C Pressure
250 preswt(ii, imat) = zero
251C Soudspeed
252 sspwt(ii, imat) = zero
253C Gruneisen coefficient
254 grun(ii, imat) = zero
255 ENDDO
256 ENDDO
257C ----------------------------------------------
258C OUTER RELAXATION LOOP
259C ----------------------------------------------
260 error(1:nel) = ep10
261C Number of elements that have not converged
262 remaining_elts = 0
263 DO ii = 1, nel
264 cv_indicator(ii) = 1
265 IF (nbmat_incell(ii) > 1) THEN
266 cv_indicator(ii) = 0
267 remaining_elts = remaining_elts + 1
268 ENDIF
269 relaxed_pressure(ii) = multi_fvm%PRES(ii + nft)
270 p0(ii) = zero
271 ENDDO
272C Loop over the fluid layers
273 DO imat = 1, nbmat
274 IF(ASSOCIATED(ix)) THEN
275 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
276 ELSE
277 local_matid = 1
278 ENDIF
279 matlaw = matparam(local_matid)%ILAW
280 pres => preswt(1:nel, imat)
281 ssp => sspwt(1:nel, imat)
282 vol => volwt(1:nel, imat)
283 rho => rhowt(1:nel, imat)
284 eint => eintwt(1:nel, imat)
285 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
286 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
287 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
288 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
289C Compute EOS
290 CALL multi_submatlaw(
291 1 1, matlaw, local_matid, nel,
292 2 eint, pres, rho, ssp,
293 3 vol, grun(1:nel,imat),pm, ipm,
294 4 npropm, npropmi, bufmat, lbuf%OFF,
295 5 lbuf%TEMP, lbuf%BFRAC, gbuf%TB, gbuf%DELTAX,
296 6 current_time, sigold, snpc, stf,
297 7 npf, tf, ebuf%VAR, nvareos,
298 7 matparam(local_matid), nvartmp_eos, ebuf%VARTMP, nummat ,
299 8 gbuf%ABURN)
300C Relaxed pressure initialization
301 DO ii = 1, nel
302 p0(ii) = p0(ii) + alphawt(ii, imat) * pres(ii)
303 ENDDO
304 ENDDO ! IMAT = 1, NBMAT
305C RHS for energy equation
306 DO imat = 1, nbmat
307 DO ii = 1, nel
308 IF (alphawt(ii, imat) /= zero) THEN
309 fe0wt(ii, imat) = speeintwt(ii, imat)
310 ENDIF
311 ENDDO
312 ENDDO
313C It seems that submaterial pressure is not initialized
314 IF (timestep == zero) THEN
315 DO ii = 1, nel
316 relaxed_pressure(ii) = p0(ii)
317 ENDDO
318 ENDIF
319
320
321C =================
322C Newton iterations
323C =================
324C Parameters
325 max_iter = 50
326 tol = em4
327 cont = .true.
328 IF (remaining_elts == 0) cont = .false.
329 iter = 0
330 DO WHILE(iter < max_iter .AND. cont)
331 iter = iter + 1
332 coef(1:nel) = zero
333 pres_incr(1:nel) = zero
334 error(1:nel) = zero
335 DO imat = 1, nbmat
336 feint => fewt(1:nel, imat)
337 fp => fpwt(1:nel, imat)
338 ssp => sspwt(1:nel, imat)
339 pres => preswt(1:nel, imat)
340 tau => tauwt(1:nel, imat)
341 rho => rhowt(1:nel, imat)
342 speeint => speeintwt(1:nel, imat)
343 DO ii = 1, nel
344 IF (cv_indicator(ii) /= 1) THEN
345 IF (alphawt(ii, imat) /= zero .AND.
346 . mfracwt(ii, imat) > em14) THEN
347 fp(ii) = pres(ii) - relaxed_pressure(ii)
348 dpde = grun(ii, imat) / tau(ii)
349 dpdtau = - (ssp(ii) - grun(ii, imat) * (pres(ii) + pshift) * tau(ii)) / tau(ii) / tau(ii)
350 coef1 = masswt(ii, imat) / (dpdtau - (relaxed_pressure(ii) + pshift) * dpde)
351 coef(ii) = coef(ii) + coef1 * (one + grun(ii, imat))
352 pres_incr(ii) = pres_incr(ii) - fp(ii) * coef1
353 ENDIF
354 ENDIF
355 ENDDO
356 ENDDO
357
358
359 !ERROR(1:NEL) = ZERO
360 DO ii = 1, nel
361 IF (cv_indicator(ii) /= 1) THEN
362 pres_incr(ii) = pres_incr(ii) / coef(ii)
363 error(ii) = abs(pres_incr(ii))
364 ENDIF
365 ENDDO
366 lim_tau(1:nel, 1:nbmat) = one
367 DO imat = 1, nbmat
368 fp => fpwt(1:nel, imat)
369 ssp => sspwt(1:nel, imat)
370 pres => preswt(1:nel, imat)
371 tau => tauwt(1:nel, imat)
372 rho => rhowt(1:nel, imat)
373 speeint => speeintwt(1:nel, imat)
374 DO ii = 1, nel
375 IF (cv_indicator(ii) /= 1) THEN
376 IF (alphawt(ii, imat) /= zero .AND.
377 . mfracwt(ii, imat) > em14) THEN
378 dpdtau = - (ssp(ii) - grun(ii, imat) * (pres(ii) + pshift) * tau(ii)) / tau(ii) / tau(ii)
379 dpde = grun(ii, imat) / tau(ii)
380 coef1 = one / (dpdtau - (relaxed_pressure(ii) + pshift) * dpde)
381 tau_incr(ii, imat) = coef1 * (pres_incr(ii) * (one + grun(ii, imat)) + fp(ii))
382 lim = one
383 IF (tau(ii) - tau_incr(ii, imat) < zero
384 . .AND. tau_incr(ii, imat) > zero) THEN
385 lim = min(one, half * tau(ii) / tau_incr(ii, imat))
386 ENDIF
387 lim_tau(ii, imat) = lim
388 ENDIF
389 ENDIF
390 ENDDO
391 ENDDO
392
393
394 DO imat = 1, nbmat
395 pres => preswt(1:nel, imat)
396 ssp => sspwt(1:nel, imat)
397 vol => volwt(1:nel, imat)
398 rho => rhowt(1:nel, imat)
399 eint => eintwt(1:nel, imat)
400 tau => tauwt(1:nel, imat)
401 DO ii = 1, nel
402 IF (cv_indicator(ii) /= 1) THEN
403 IF (alphawt(ii, imat) /= zero .AND.
404 . mfracwt(ii, imat) > em14) THEN
405 lim = minval(lim_tau(ii, 1:nbmat))
406 tau(ii) = tau(ii) - lim * tau_incr(ii, imat)
407
408 rho(ii) = one / tau(ii)
409 alphawt(ii, imat) = masswt(ii, imat) * tau(ii)
410 vol(ii) = alphawt(ii, imat) * gbuf%VOL(ii)
411 ENDIF
412 ENDIF
413 ENDDO
414 ENDDO
415 DO ii = 1, nel
416 IF (cv_indicator(ii) /= 1) THEN
417 relaxed_pressure(ii) = relaxed_pressure(ii) - pres_incr(ii)
418 ENDIF
419 ENDDO
420 DO imat = 1, nbmat
421 eint => eintwt(1:nel, imat)
422 rho => rhowt(1:nel, imat)
423 speeint => speeintwt(1:nel, imat)
424 tau => tauwt(1:nel, imat)
425 DO ii = 1, nel
426 IF (cv_indicator(ii) /= 1) THEN
427 IF (alphawt(ii, imat) /= zero .AND.
428 . mfracwt(ii, imat) > em14) THEN
429C IF (ALPHAWT(II, IMAT) /= ZERO .AND. GRUN(II, IMAT) /= ZERO) THEN
430 speeint(ii) = fe0wt(ii, imat) -
431 . (relaxed_pressure(ii) + pshift) * (tau(ii) - tau0wt(ii, imat))
432 eint(ii) = rho(ii) * speeint(ii)
433 ENDIF
434 ENDIF
435 ENDDO
436 ENDDO
437 DO ii = 1, nel
438 IF (cv_indicator(ii) /= 1) THEN
439 IF (error(ii) < tol * (one + abs(relaxed_pressure(ii)))) THEN
440 cv_indicator(ii) = 1
441 remaining_elts = remaining_elts - 1
442 ENDIF
443 ENDIF
444 ENDDO
445 IF (remaining_elts == 0) THEN
446 cont = .false.
447 ENDIF
448C Equation of state call
449 DO imat = 1, nbmat
450 IF(ASSOCIATED(ix)) THEN
451 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
452 ELSE
453 local_matid = 1
454 ENDIF
455 matlaw = matparam(local_matid)%ILAW
456 pres => preswt(1:nel, imat)
457 ssp => sspwt(1:nel, imat)
458 vol => volwt(1:nel, imat)
459 rho => rhowt(1:nel, imat)
460 eint => eintwt(1:nel, imat)
461 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
462 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
463 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
464 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
465C Compute EOS
466 CALL multi_submatlaw(
467 1 0, matlaw, local_matid, nel,
468 2 eint, pres, rho, ssp,
469 3 vol, grun(1:nel,imat),pm, ipm,
470 4 npropm, npropmi, bufmat, lbuf%OFF,
471 5 lbuf%TEMP, lbuf%BFRAC, gbuf%TB, gbuf%DELTAX,
472 6 current_time, sigold, snpc, stf,
473 7 npf, tf, ebuf%VAR, nvareos,
474 8 matparam(local_matid), nvartmp_eos, ebuf%VARTMP , nummat ,
475 9 gbuf%ABURN)
476 ENDDO ! IMAT = 1, NBMAT
477
478 ENDDO !ITER = 1, MAX_ITER
479 DO ii = 1, nel
480 DO imat = 1, nbmat
481 IF (eintwt(ii, imat) < zero) THEN
482 !!!PRINT*, "OUPS1"
483 eintwt(ii, imat) = -eintwt(ii, imat)
484 ENDIF
485 ENDDO
486 ENDDO
487 IF (remaining_elts /= 0 .AND. idbg) THEN
488 print*, "*** CV PB IN PRESSURE EQUILIBRIUM", " CYCLE", ncycle, "GROUPE", ng
489 DO ii = 1, nel
490 IF (cv_indicator(ii) == 0 .AND. ASSOCIATED(ix)) THEN
491 WRITE(*, '(A,I10,A,I10)') "LOCAL_ID ", ii + nft, " GLOBAL_ID ", ix(nix, ii)
492 DO imat = 1, nbmat
493 WRITE(*,*) "MAT ", imat, ", alph : ", alphawt(ii, imat), ", PRES : ", preswt(ii, imat)
494 . , ", EINT : ", eintwt(ii, imat)
495 ENDDO
496 CONTINUE
497 ENDIF
498 ENDDO
499 ENDIF
500C ----------------------------------------------
501C END OF OUTER RELAXATION LOOP
502C ----------------------------------------------
503
504C Loop over the fluid layers
505 DO imat = 1, nbmat
506 IF(ASSOCIATED(ix)) THEN
507 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
508 ELSE
509 local_matid = 1
510 ENDIF
511 matlaw = matparam(local_matid)%ILAW
512 pres => preswt(1:nel, imat)
513 ssp => sspwt(1:nel, imat)
514 vol => volwt(1:nel, imat)
515 rho => rhowt(1:nel, imat)
516 eint => eintwt(1:nel, imat)
517
518 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
519 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
520 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
521
522 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
523 multi_fvm%PHASE_EINT(imat, 1 + nft : nel + nft) = eintwt(1:nel, imat)
524 multi_fvm%PHASE_RHO(imat, 1 + nft : nel + nft) = rhowt(1:nel, imat)
525 DO ii = 1, nel
526 i = ii + nft
527 multi_fvm%PHASE_ALPHA(imat, i) = volwt(ii, imat) / gbuf%VOL(ii)
528 ENDDO
529C Compute EOS
530 CALL multi_submatlaw(
531 1 0, matlaw, local_matid, nel,
532 2 eint, pres, rho, ssp,
533 3 vol, grun(1:nel,imat),pm, ipm,
534 4 npropm, npropmi, bufmat, lbuf%OFF,
535 5 lbuf%TEMP, lbuf%BFRAC, gbuf%TB, gbuf%DELTAX,
536 6 current_time, sigold, snpc, stf,
537 7 npf, tf, ebuf%VAR, nvareos,
538 8 matparam(local_matid), nvartmp_eos, ebuf%VARTMP, nummat ,
539 9 gbuf%ABURN)
540
541 ENDDO ! IMAT = 1, NBMAT
542
543C Equilibrium densities have been found,
544C ----------------------------------------------
545C IT IS NECESSARY NOW TO RECOVER THE GLOBAL ENERGY
546C THAT IS ENSURING THAT E = E1 + E2 + ...
547C ----------------------------------------------
548C We're solving :
549C SUM ALPHA_I EINT_I - EINT = 0
550C P1 - P = 0
551C ...
552C PN - P = 0
553C With unknowns P, EINT1, EINT2, ..., EINT_N
554C Initial guess for P is the relaxed pressure coming from the
555C Previous iterative system
556C Initial guesses for submaterial internal energy also come
557C from the previous system
558 error(1:nel) = ep10
559C Number of elements that have not converged
560 remaining_elts = 0
561 DO ii = 1, nel
562 i = ii + nft
563 cv_indicator(ii) = 1
564 one_over_grun = one
565 DO imat = 1, nbmat
566 IF (multi_fvm%PHASE_ALPHA(imat, i) > zero .AND.
567 . mfracwt(ii, imat) > em14) THEN
568 one_over_grun = one_over_grun * grun(ii, imat)
569 ENDIF
570 ENDDO
571 IF (nbmat_incell(ii) > 1 .AND. one_over_grun > zero) THEN
572C IF (NBMAT_INCELL(II) > 1) THEN
573 cv_indicator(ii) = 0
574 remaining_elts = remaining_elts + 1
575 ENDIF
576 ENDDO
577
578 cont = .true.
579 IF (remaining_elts == 0) cont = .false.
580 iter = 0
581 tole = em06
582 max_iter = 10
583 DO WHILE(iter < max_iter .AND. cont)
584 iter = iter + 1
585C Compute new pressure
586 ! Store coefficients
587 coef(1:nel) = zero
588 ! Function to cancel
589 fe(1:nel) = - multi_fvm%EINT(1 + nft : nel + nft)
590 pres_incr(1:nel) = zero
591 DO imat = 1, nbmat
592 pres => preswt(1:nel, imat)
593 vol => volwt(1:nel, imat)
594 eint => eintwt(1:nel, imat)
595 fp => fpwt(1:nel, imat)
596 DO ii = 1, nel
597 IF (vol(ii) > zero .AND.
598 . mfracwt(ii, imat) > em14) THEN
599 IF (cv_indicator(ii) /= 1) THEN
600 one_over_grun = zero
601 IF (grun(ii, imat) > zero) THEN
602 one_over_grun = one / grun(ii, imat)
603 ENDIF
604 vol_frac = alphawt(ii, imat)
605
606 fp(ii) = pres(ii) - relaxed_pressure(ii)
607 IF (one_over_grun > zero) THEN
608 fe(ii) = fe(ii) + vol_frac * eint(ii)
609 ENDIF
610
611 pres_incr(ii) = pres_incr(ii) + vol_frac * one_over_grun * fp(ii)
612 coef(ii) = coef(ii) + vol_frac * one_over_grun
613
614 ENDIF
615 ENDIF
616 ENDDO ! II = 1, NEL
617 ENDDO ! IMAT = 1, NBMAT
618 DO ii = 1, nel
619 IF (cv_indicator(ii) /= 1) THEN
620 IF (coef(ii) /= zero) THEN
621! Pressure increment
622 pres_incr(ii) = (fe(ii) - pres_incr(ii)) / coef(ii)
623 relaxed_pressure(ii) = relaxed_pressure(ii) - pres_incr(ii)
624! Compute ERROR
625 error(ii) = abs(pres_incr(ii))
626 ENDIF
627 ENDIF
628 ENDDO
629C Compute New energies
630 DO imat = 1, nbmat
631 pres => preswt(1:nel, imat)
632 vol => volwt(1:nel, imat)
633 rho => rhowt(1:nel, imat)
634 eint => eintwt(1:nel, imat)
635 fp => fpwt(1:nel, imat)
636 DO ii = 1, nel
637 IF (cv_indicator(ii) /= 1) THEN
638 IF (vol(ii) > zero .AND. grun(ii, imat) > zero .AND.
639 . mfracwt(ii, imat) > em14) THEN
640 lim = one
641 IF ((fp(ii) + pres_incr(ii)) > zero) THEN
642 lim = max(zero, min (one, half * eint(ii) * grun(ii, imat)
643 . / (fp(ii) + pres_incr(ii))))
644 ENDIF
645 eint(ii) = eint(ii) - lim * (fp(ii) + pres_incr(ii)) / grun(ii, imat)
646 ENDIF
647 ENDIF
648 ENDDO ! II = 1, NEL
649C Call equation of state for the submaterial
650 IF(ASSOCIATED(ix)) THEN
651 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
652 ELSE
653 local_matid = 1
654 ENDIF
655 matlaw = matparam(local_matid)%ILAW
656 ssp => sspwt(1:nel, imat)
657 vol => volwt(1:nel, imat)
658 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
659 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
660 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
661 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
662 CALL multi_submatlaw(
663 1 0, matlaw, local_matid, nel,
664 2 eint, pres, rho, ssp,
665 3 vol, grun(1:nel,imat),pm, ipm,
666 4 npropm, npropmi, bufmat, lbuf%OFF,
667 5 lbuf%TEMP, lbuf%BFRAC, gbuf%TB, gbuf%DELTAX,
668 6 current_time, sigold, snpc, stf,
669 7 npf, tf, ebuf%VAR, nvareos,
670 8 matparam(local_matid),nvartmp_eos, ebuf%VARTMP, nummat ,
671 9 gbuf%ABURN)
672 ENDDO ! IMAT = 1, NBMAT
673 DO ii = 1, nel
674 IF (cv_indicator(ii) /= 1) THEN
675 IF (error(ii) < tol * (one + abs(relaxed_pressure(ii)))) THEN
676 cv_indicator(ii) = 1
677 remaining_elts = remaining_elts - 1
678 ENDIF
679 ENDIF
680 ENDDO
681 IF (remaining_elts == 0) THEN
682 cont = .false.
683 ENDIF
684 ENDDO ! WHILE(ITER < MAX_ITER .AND. CONT)
685 DO ii = 1, nel
686 DO imat = 1, nbmat
687 IF (eintwt(ii, imat) < zero) THEN
688C PRINT*, "OUPS2"
689 ENDIF
690 ENDDO
691 ENDDO
692 IF (remaining_elts /= 0 .AND. cont .AND. idbg) THEN
693 print*, "*** CV PB IN ENERGY RESET", "CYCLE", ncycle, "GROUPE", ng
694 ENDIF
695C ----------------------------------------------
696C END OF ENERGY RESET
697C THAT IS ENSURING THAT E = E1 + E2 + ...
698C ----------------------------------------------
699 gbuf%SIG(1:nel) = zero
700C Store final values
701 multi_fvm%SOUND_SPEED(1 + nft:nel + nft) = zero
702 DO imat = 1, nbmat
703 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
704 multi_fvm%PHASE_EINT(imat, 1 + nft : nel + nft) = eintwt(1:nel, imat)
705 multi_fvm%PHASE_RHO(imat, 1 + nft : nel + nft) = rhowt(1:nel, imat)
706C Pressure
707 multi_fvm%PHASE_PRES(imat, 1 + nft : nel + nft) = preswt(1:nel, imat)
708 IF(ASSOCIATED(ix)) THEN
709 local_matid = matparam(ix(1,1))%MULTIMAT%MID(imat)
710 ELSE
711 local_matid = 1
712 ENDIF
713
714 matlaw = matparam(local_matid)%ILAW
715C Global Pressure: volumic fraction ponderation
716 DO ii = 1, nel
717 i = ii + nft
718 vol_frac = multi_fvm%PHASE_ALPHA(imat, i)
719 gbuf%SIG(ii) = gbuf%SIG(ii) - preswt(ii, imat) * vol_frac
720 multi_fvm%SOUND_SPEED(ii + nft) = multi_fvm%SOUND_SPEED(ii + nft) +
721 . vol_frac * rhowt(ii, imat) * max(sspwt(ii, imat), em20)
722 IF (matlaw == 5) THEN
723 multi_fvm%BFRAC(imat, i) = lbuf%BFRAC(ii)
724 ELSE
725 multi_fvm%BFRAC(imat, i) = zero
726 ENDIF
727 ENDDO
728 ENDDO
729c$$$C GLobal density (should not have changed except for numerical errors
730c$$$ DO II = 1, NEL
731c$$$ I = II + NFT
732c$$$ MULTI_FVM%RHO(I) =ZERO
733c$$$ DO IMAT = 1, NBMAT
734c$$$ MULTI_FVM%RHO(I) = MULTI_FVM%RHO(I) +
735c$$$ . MULTI_FVM%PHASE_ALPHA(IMAT, I) * MULTI_FVM%PHASE_RHO(IMAT, I)
736c$$$ ENDDO
737c$$$ ENDDO
738 gbuf%SIG(1 + 1 * nel : 2 * nel) = gbuf%SIG(1 + 0 * nel : 1 * nel)
739 gbuf%SIG(1 + 2 * nel : 3 * nel) = gbuf%SIG(1 + 0 * nel : 1 * nel)
740C Global Sound speed
741 DO ii = 1, nel
742 i = ii + nft
743 multi_fvm%SOUND_SPEED(ii + nft) = sqrt(multi_fvm%SOUND_SPEED(ii + nft)
744 . / multi_fvm%RHO(i))
745 multi_fvm%PRES(ii + nft) = -gbuf%SIG(ii)
746 ENDDO
747C Global Temperature
748 !we must introduce an iterative solver : solve "find T such as e_global+P/rho = int(Cp_global(T) dT)"
749 !this requires also temperature solving from all EoS (improvement to plan)
750 gbuf%TEMP(1:nel) = zero
751 DO imat=1,nbmat
752 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
753 DO ii = 1, nel
754 gbuf%TEMP(ii) = gbuf%TEMP(ii) + mfracwt(ii, imat)*lbuf%TEMP(ii)
755 ENDDO
756 ENDDO
757 ENDIF ! NBMAT == 1
758 ENDIF ! MTN == 151
759 ENDDO ! ng = itask + 1, ngroup, nthread
760C-----------------------------------------------
761C E n d o f s u b r o u t i n e
762C-----------------------------------------------
763 END SUBROUTINE multi_pressure_equilibrium
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_pressure_equilibrium
subroutine multi_pressure_equilibrium(timestep, elbuf_tab, iparg, itask, ixs, ixq, ixtg, pm, ipm, multi_fvm, current_time, bufmat, npf, tf, nummat, matparam)
Definition multi_pressure_equilibrium.F:38
initbuf_mod
Definition initbuf.F:160
multi_submatlaw_mod
Definition multi_submatlaw.F:34
multi_submatlaw_mod::multi_submatlaw
subroutine multi_submatlaw(iflag, matlaw, local_matid, nel, eint, pres, rho, ssp, vol, grun, pm, ipm, npropm, npropmi, bufmat, off, theta, burnfrac, burntime, deltax, current_time, sigold, snpc, stf, npf, tf, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos, nummat, aburn)
Definition multi_submatlaw.F:63