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