OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
inigrav_eos.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!|| inigrav_eos ../starter/source/initial_conditions/inigrav/inigrav_eos.F
25!||--- called by ------------------------------------------------------
26!|| inigrav_load ../starter/source/initial_conditions/inigrav/inigrav_load.F
27!||--- calls -----------------------------------------------------
28!||--- uses -----------------------------------------------------
29!||====================================================================
30 SUBROUTINE inigrav_eos(NELG, NEL, NG, MATID, IPM, GRAV0, RHO0, DEPTH, PM, BUFMAT,
31 . ELBUF_TAB,PSURF,LIST,PRES,IMAT,MLW,NPF,TF,NUMMAT,MAT_PARAM)
32C-----------------------------------------------
33C M o d u l e s
34C-----------------------------------------------
35 USE elbufdef_mod
36 USE matparam_def_mod
37 USE eosmain_mod , ONLY : eosmain
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C-----------------------------------------------
43C C o m m o n B l o c k s
44C-----------------------------------------------
45#include "param_c.inc"
46#include "com01_c.inc"
47#include "tabsiz_c.inc"
48C-----------------------------------------------
49C D u m m y A r g u m e n t s
50C-----------------------------------------------
51 INTEGER, INTENT(IN) :: NUMMAT
52 INTEGER, INTENT(IN) :: NEL, NG, MATID, IPM(NPROPMI, *), LIST(NEL),NELG, IMAT, MLW
53 my_real, INTENT(IN) :: grav0, depth(*), pm(npropm, *), bufmat(*),psurf
54 my_real, INTENT(INOUT) :: pres(nel)
55 TYPE(elbuf_struct_), DIMENSION(NGROUP), TARGET, INTENT(IN) :: ELBUF_TAB
56 my_real, INTENT(IN) :: rho0
57 INTEGER,INTENT(IN)::NPF(SNPC)
58 my_real,INTENT(IN)::tf(stf)
59 TYPE(matparam_struct_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
60C-----------------------------------------------
61C L o c a l V a r i a b l e s
62C-----------------------------------------------
63 INTEGER :: I,K
64 INTEGER :: EOSTYPE, ITER, MAX_ITER, MAT(NEL), REMAINING_ELTS, NBMAT, NVAREOS, NVARTMP_EOS
65 my_real :: bid(nel) !bfrac not used int his context
66 my_real :: tol, incr
67 my_real :: func(nel), dfunc(nel), error(nel), rhozero(nel), rho(nel),
68 . mu(nel), mu2(nel), espe(nel), dvol(nel), df(nel), vol(nel), psh(nel),
69 . dpdrho(nel), dpdm(nel), dpde(nel), theta(nel), ecold(nel), p0(nel),
70 . grun(nel),pref(nel), off(nel), eint(nel), sig(nel,6),rho_old,v0(nel)
71
72 LOGICAL :: CONT, CONVERGED(NEL)
73 TYPE(g_bufel_), POINTER :: GBUF
74 TYPE(l_bufel_), POINTER :: LBUF
75 TYPE(buf_eos_), POINTER :: EBUF
76C-----------------------------------------------
77C S o u r c e L i n e s
78C-----------------------------------------------
79
80C LIST IS SUBGROUP TO TREAT : ONLY ELEM WITH RELEVANT PARTS ARE KEPT
81C NEL IS ISEZ OF LIST
82C NELG IS SIZE OF ORIGINAL GROUP : needed to shift indexes in GBUF%SIG & MBUF%VAR
83C IMAT : relevant for multifluid eos. imat=0 => monomaterial law
84
85 gbuf => elbuf_tab(ng)%GBUF
86 NULLIFY(lbuf)
87 NULLIFY(ebuf)
88 ebuf => elbuf_tab(ng)%BUFLY(1)%EOS(1, 1, 1)
89 nvareos = elbuf_tab(ng)%BUFLY(1)%NVAR_EOS
90 nvartmp_eos = elbuf_tab(ng)%BUFLY(1)%NVARTMP_EOS
91 IF(imat>0)THEN
92 lbuf => elbuf_tab(ng)%BUFLY(imat)%LBUF(1, 1, 1)
93 ebuf => elbuf_tab(ng)%BUFLY(imat)%EOS(1, 1, 1)
94 nvareos = elbuf_tab(ng)%BUFLY(imat)%NVAR_EOS
95 nvartmp_eos = elbuf_tab(ng)%BUFLY(imat)%NVARTMP_EOS
96 ENDIF
97
98 eostype = mat_param(matid)%IEOS
99 eostype = ipm(4, matid) ! still used temporarily in this section
100 iter = 0
101 max_iter = 30
102 cont =.true.
103 error(1:nel) = ep30
104 tol = em4
105 converged(1:nel) =.false.
106
107 !===================================================================!
108 ! BUFFER INITIALIZATION !
109 !===================================================================!
110 ! local array instead of pointers since they may be reshaped with a
111 ! subgroup (1:NEL) c (1:NELG) depending on element PART (possibly several PARTS in a given group)
112 DO k = 1, nel
113 rhozero(k) = pm(1,matid)
114 psh(k) = pm(88, matid)
115 mat(k) = matid
116 pref(k) = psurf - rho0 * grav0 * depth(k) !RHO in : RHO0_mixture (gravity model applied to the mixture not to the single submaterial)
117 vol(k) = one ! Real volume value not useful here, just an indicator : ONE = compute, ZERO = do not compute
118 mu(k) = zero
119 sig(k,1) = -pm(31,matid)
120 sig(k,2) = -pm(31,matid)
121 sig(k,3) = -pm(31,matid)
122 sig(k,4) = zero
123 sig(k,5) = zero
124 sig(k,6) = zero
125 ENDDO
126
127 IF(imat==0)THEN ! WORK IN GBUF
128 DO k = 1, nel
129 i = list(k)
130 rho(k) = gbuf%RHO(k)
131 off(k) = gbuf%OFF(i)
132 eint(k) = gbuf%EINT(i)
133 v0(k) = gbuf%VOL(i)
134 ENDDO
135 ELSE ! THEN IMAT>0 WORK IN LBUF
136 DO k = 1, nel
137 i = list(k)
138 rho(k) = pm(1,matid) !submaterial
139 off(k) = lbuf%OFF(i)
140 eint(k) = lbuf%EINT(i)
141 v0(k) = lbuf%VOL(i)
142 ENDDO
143 ENDIF
144 !EINT=Eint/V in this context with law 151
145 !===================================================================!
146 ! iterations !
147 !===================================================================!
148 remaining_elts = nel !Number of elements for which convergence is not achieved
149 DO WHILE(cont .AND. iter < max_iter) !Start nonlinear solver loop
150 iter = iter + 1
151 DO k = 1, nel
152 mu(k) = rho(k) / rhozero(k) - one
153 df(k) = rhozero(k) / rho(k)
154 eint(k) = eint(k) / df(k) !with IFLAG==2 we need ESPE=Df*EINT
155 ENDDO
156 CALL eosmain(2 , nel , eostype , pm , off , eint,
157 . rho , rhozero , mu , mu2 , espe ,
158 . dvol , df , vol , mat , psh ,
159 . pres , dpdm , dpde , theta ,
160 . bufmat, sig , mu , mlw ,
161 . npf , tf , ebuf%VAR , nvareos, mat_param(matid),
162 . bid , nvartmp_eos, ebuf%VARTMP)
163 DO k=1,nel
164 eint(k)=eint(k)*df(k) ! set back to Eint/V
165 ENDDO
166 !===================================================================!
167 ! CHECK CONVERGENCY !
168 !===================================================================!
169 !loop over retained element in the current group I \in LIST(1:NEL)
170 DO k = 1, nel
171 IF (.NOT. converged(k)) THEN ! No need to take the square root, only the squared sound speed is needed here
172 dpdrho(k) = dpdm(k) / rhozero(k) !DPDM is total derivative (variation along isentropic path)
173 !GRUN(K) = DPDE(K) / (ONE + MU(K))
174 func(k) = pres(k)-pref(k)
175 dfunc(k) = dpdrho(k) !- GRUN(K) * PRES(K) / RHO(K)
176 incr = - func(k) / dfunc(k)
177 rho_old = rho(k)
178 rho(k) = rho(k) + incr
179 error(k)= max( abs(incr)/rhozero(k) , abs(func(k)) )
180
181 incr = - (pres(k)+psh(k)+half*dpdrho(k)*incr)*(rhozero(k)/rho(k)-rhozero(k)/rho_old)
182 eint(k) = eint(k) + incr !EINT is in fact Eint/V with law 151 here
183 IF(dpde(k)==zero)eint(k)=zero
184
185 IF (error(k) < tol .OR. iter == max_iter) THEN
186 converged(k) = .true.
187 remaining_elts = remaining_elts - 1
188 vol(k) = zero
189 i = list(k)
190 df(k) = rhozero(i) / rho(k)
191 IF(imat==0)THEN
192 gbuf%RHO(i) = rho(k)
193 gbuf%SIG(i ) = - pres(k)
194 gbuf%SIG(i + nelg) = - pres(k)
195 gbuf%SIG(i + 2 * nelg) = - pres(k)
196 gbuf%EINT(i) = eint(k)/df(k) !Eint/V with law 151 here => EINT*DF == Eint/V0
197 ELSE
198 lbuf%RHO(i) = rho(k)
199 lbuf%SIG(i ) = - pres(k)
200 lbuf%SIG(i + nelg) = - pres(k)
201 lbuf%SIG(i + 2 * nelg) = - pres(k)
202 lbuf%EINT(i) = eint(k)/df(k) !Eint/V with law 151 here
203 ENDIF
204 ENDIF
205 ENDIF
206 ENDDO
207 cont = (remaining_elts /= 0)
208 ENDDO ! WHILE
209
210
211 END SUBROUTINE inigrav_eos
my_real
#define my_real
Definition cppsort.cpp:32
inigrav_eos
subroutine inigrav_eos(nelg, nel, ng, matid, ipm, grav0, rho0, depth, pm, bufmat, elbuf_tab, psurf, list, pres, imat, mlw, npf, tf, nummat, mat_param)
Definition inigrav_eos.F:32
max
#define max(a, b)
Definition macros.h:21
eosmain_mod
Definition eosmain.F:33
eosmain_mod::eosmain
subroutine eosmain(iflag, nel, eostyp, pm, off, eint, rho, rho0, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, theta, bufmat, sig, mu_bak, mlw, npf, tf, vareos, nvareos, mat_param, bfrac, nvartmp, vartmp)
Definition eosmain.F:80