OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
inigrav_m51.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_m51 ../starter/source/initial_conditions/inigrav/inigrav_m51.f
25!||--- called by ------------------------------------------------------
26!|| inigrav_load ../starter/source/initial_conditions/inigrav/inigrav_load.F
27!||--- calls -----------------------------------------------------
28!|| polysolve ../starter/source/initial_conditions/inigrav/inigrav_m51.F
29!||--- uses -----------------------------------------------------
30!||====================================================================
31 SUBROUTINE inigrav_m51(NELG, NEL , NG , MATID, IPM, GRAV0, DEPTH, PM, BUFMAT, ELBUF_TAB,
32 . PSURF,LIST, ALE_CONNECTIVITY, NV46 , IX , NIX , NFT , BUFMATG, IPARG)
33C-----------------------------------------------
34C M o d u l e s
35C-----------------------------------------------
36 USE elbufdef_mod
37 USE ale_connectivity_mod
38 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C-----------------------------------------------
44C C o m m o n B l o c k s
45C-----------------------------------------------
46! NPROPMI, NPROPM
47#include "param_c.inc"
48! NGROUP
49#include "com01_c.inc"
50C-----------------------------------------------
51C D u m m y A r g u m e n t s
52C-----------------------------------------------
53 INTEGER, INTENT(IN) :: NEL, NG, MATID, IPM(NPROPMI, *),LIST(NEL),NELG,IX(NIX,*),NFT,NIX,
54 . iparg(nparg,ngroup)
55 my_real, INTENT(IN) :: grav0, depth(*), pm(npropm, *), bufmat(*),psurf,bufmatg(*)
56 TYPE(elbuf_struct_), DIMENSION(NGROUP), TARGET, INTENT(IN) :: ELBUF_TAB
57 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
58C-----------------------------------------------
59C L o c a l V a r i a b l e s
60C-----------------------------------------------
61 INTEGER :: I, K1, K2, K3, K4, NUVAR,K,J,IFORM,IFORMv,NV46,IV,IADBUF,ML,NGv,KTY,KLT,N,MFT,IS,NELGv
62 my_real :: pgrav, rho, gam, rho0,
63 . alpha1, alpha2, alpha3, alpha4,
64 . c01, c11, c21, c31, c41, c51,
65 . c02, c12, c22, c32, c42, c52,
66 . c03, c13, c23, c33, c43, c53,
67 . c04, pext,
68 . rho10, rho20, rho30, rho40, rho1, rho2, rho3, rho4, mu,
69 . eint1, eint2, eint3, eint4, vol, vol1, vol2, vol3,
70 . eint10, eint20, eint30, eint40
71 TYPE(G_BUFEL_), POINTER :: GBUF
72 TYPE(buf_mat_) ,POINTER :: MBUF , MBUFv
73 INTEGER :: IAD
74C-----------------------------------------------
75C S o u r c e L i n e s
76C-----------------------------------------------
77
78C LIST IS SUBGROUP TO TREAT : ONLY ELEM WITH RELEVANT PARTS ARE KEPT
79C NEL IS ISEZ OF LIST
80C NELG IS SIZE OF ORIGINAL GROUP : needed to shift indexes in GBUF%SIG & MBUF%VAR
81
82 !Global buffer
83 gbuf => elbuf_tab(ng)%GBUF
84 !Material buffer
85 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
86
87 iform = nint(bufmat(31))
88 k1 = m51_n0phas + (1 - 1) * m51_nvphas
89 k2 = m51_n0phas + (2 - 1) * m51_nvphas
90 k3 = m51_n0phas + (3 - 1) * m51_nvphas
91 k4 = m51_n0phas + (4 - 1) * m51_nvphas
92 nuvar = ipm(8, matid)
93 pext = bufmat(8)
94 IF(iform /= 6)THEN
95 !all value are in UPARAM=>BUFMAT()
96 !Initial densities
97 rho10 = bufmat(9)
98 rho20 = bufmat(10)
99 rho30 = bufmat(11)
100 rho40 = bufmat(12)
101 !Initial energies
102 eint10 = bufmat(32)
103 eint20 = bufmat(33)
104 eint30 = bufmat(34)
105 eint40 = bufmat(48)
106 !Internal energies
107 ! EINT1 = MBUF%VAR(I + (K1 + 8 - 1) * NELG)
108 ! EINT2 = MBUF%VAR(I + (K2 + 8 - 1) * NELG)
109 ! EINT3 = MBUF%VAR(I + (K3 + 8 - 1) * NELG)
110 ! EINT4 = MBUF%VAR(I + (K4 + 8 - 1) * NELG)
111 !Material 1
112 c01 = bufmat(35)
113 c11 = bufmat(12)
114 c21 = bufmat(15)
115 c31 = bufmat(18)
116 c41 = bufmat(22)
117 c51 = bufmat(25)
118 !Material 2
119 c02 = bufmat(36)
120 c12 = bufmat(13)
121 c22 = bufmat(16)
122 c32 = bufmat(20)
123 c42 = bufmat(23)
124 c52 = bufmat(26)
125 !Material 3
126 c03 = bufmat(37)
127 c13 = bufmat(14)
128 c23 = bufmat(17)
129 c33 = bufmat(21)
130 c43 = bufmat(24)
131 c53 = bufmat(27)
132 !Material 4
133 c04 = bufmat(49)
134 ENDIF
135 DO k = 1, nel
136 i = list(k)
137 IF(iform==6)THEN
138 !value are not in buffer UPARAM=>BUFMAT but are cell dependent
139 ml = 0
140 iformv = 0
141 iad = ale_connectivity%ee_connect%iad_connect(i + nft)
142 DO j=1,ale_connectivity%ee_connect%iad_connect(i + nft + 1) - ale_connectivity%ee_connect%iad_connect(i + nft)
143 iv = ale_connectivity%ee_connect%connected(iad + j - 1)
144 iformv = 1000
145 IF(iv/=0) ml = nint(pm(19,ix(1,iv))) !ix is here local I and not IV
146 IF(iv/=0) iadbuf = ipm(7,ix(1,iv))
147 IF(iv/=0.AND.ml==51) iformv = bufmatg(iadbuf+31-1)
148 IF(ml==51.AND.iformv<=1)EXIT
149 ENDDO
150 DO n=1,ngroup
151 kty = iparg(5,n)
152 klt = iparg(2,n)
153 mft = iparg(3,n)
154 IF (kty==1 .AND. iv<=klt+mft) EXIT
155 ENDDO
156 IF (kty/=1 .OR. iv>klt+mft) cycle
157 ngv = n
158 mbufv => elbuf_tab(ngv)%BUFLY(1)%MAT(1,1,1)
159 nelgv = klt
160 is = iv-mft
161 rho10 = bufmatg(iadbuf-1+09)
162 rho20 = bufmatg(iadbuf-1+10)
163 rho30 = bufmatg(iadbuf-1+11)
164 rho40 = bufmatg(iadbuf-1+12)
165 !energies
166 eint1 = bufmatg(iadbuf-1+32)
167 eint2 = bufmatg(iadbuf-1+33)
168 eint3 = bufmatg(iadbuf-1+34)
169 eint4 = bufmatg(iadbuf-1+48)
170 !Material 1
171 c01 = bufmatg(iadbuf-1+35)
172 c11 = bufmatg(iadbuf-1+12)
173 c21 = bufmatg(iadbuf-1+15)
174 c31 = bufmatg(iadbuf-1+18)
175 c41 = bufmatg(iadbuf-1+22)
176 c51 = bufmatg(iadbuf-1+25)
177 !Material 2
178 c02 = bufmatg(iadbuf-1+36)
179 c12 = bufmatg(iadbuf-1+13)
180 c22 = bufmatg(iadbuf-1+16)
181 c32 = bufmatg(iadbuf-1+20)
182 c42 = bufmatg(iadbuf-1+23)
183 c52 = bufmatg(iadbuf-1+26)
184 !Material 3
185 c03 = bufmatg(iadbuf-1+37)
186 c13 = bufmatg(iadbuf-1+14)
187 c23 = bufmatg(iadbuf-1+17)
188 c33 = bufmatg(iadbuf-1+21)
189 c43 = bufmatg(iadbuf-1+24)
190 c53 = bufmatg(iadbuf-1+27)
191 !Material 4
192 c04 = bufmatg(iadbuf-1+49)
193 !vol frac
194 alpha1 = mbufv%VAR(is + (k1 + 23 - 1) * nelgv)
195 alpha2 = mbufv%VAR(is + (k2 + 23 - 1) * nelgv)
196 alpha3 = mbufv%VAR(is + (k3 + 23 - 1) * nelgv)
197 alpha4 = mbufv%VAR(is + (k4 + 23 - 1) * nelgv)
198 ELSE
199 eint1 = eint10
200 eint2 = eint20
201 eint3 = eint30
202 eint4 = eint40
203 !Volumic fractions
204 alpha1 = mbuf%VAR(i + (k1 + 23 - 1) * nelg)
205 alpha2 = mbuf%VAR(i + (k2 + 23 - 1) * nelg)
206 alpha3 = mbuf%VAR(i + (k3 + 23 - 1) * nelg)
207 alpha4 = mbuf%VAR(i + (k4 + 23 - 1) * nelg)
208 ENDIF
209
210 !Mean initial density
211 rho0 = alpha1 * rho10 + alpha2 * rho20 + alpha3 * rho30 + alpha4 * rho40
212 !Hydrostatic pressure
213 pgrav = psurf - rho0 * grav0 * depth(k)
214 !Solve for partial densities
215 vol = gbuf%VOL(i)
216 vol1 = alpha1*vol
217 vol2 = alpha2*vol
218 vol3 = alpha3*vol
219C
220 CALL polysolve(c01, c11, c21, c31, c41, c51, pgrav, eint1, mu, rho10, vol1,pext,ix(11,i+nft))
221 rho1 = rho10 * (mu + one)
222 CALL polysolve(c02, c12, c22, c32, c42, c52, pgrav, eint2, mu, rho20, vol2,pext,ix(11,i+nft))
223 rho2 = rho20 * (mu + one)
224 CALL polysolve(c03, c13, c23, c33, c43, c53, pgrav, eint3, mu, rho30, vol3,pext,ix(11,i+nft))
225 rho3 = rho30 * (mu + one)
226 !Store partial densities
227 mbuf%VAR(i + (k1 + 09 - 1) * nelg) = rho1 !initialize rho(t=0)
228 mbuf%VAR(i + (k2 + 09 - 1) * nelg) = rho2
229 mbuf%VAR(i + (k3 + 09 - 1) * nelg) = rho3
230 mbuf%VAR(i + (k1 + 12 - 1) * nelg) = rho1 !initialize rhoOLD(t=0)
231 mbuf%VAR(i + (k2 + 12 - 1) * nelg) = rho2
232 mbuf%VAR(i + (k3 + 12 - 1) * nelg) = rho3
233 mbuf%VAR(i + (k1 + 20 - 1) * nelg) = rho1 !keep rho0 (initial state is element dependent)
234 mbuf%VAR(i + (k2 + 20 - 1) * nelg) = rho2
235 mbuf%VAR(i + (k3 + 20 - 1) * nelg) = rho3
236 !Store pressures
237 mbuf%VAR(i + (k1 + 18 - 1) * nelg) = pgrav !P(t=0)
238 mbuf%VAR(i + (k2 + 18 - 1) * nelg) = pgrav
239 mbuf%VAR(i + (k3 + 18 - 1) * nelg) = pgrav
240 mbuf%VAR(i + (k4 + 18 - 1) * nelg) = pgrav
241 !Internal energies
242 IF (rho10 /= zero) THEN
243 mbuf%VAR(i + (k1 + 08 - 1) * nelg) = eint1 * rho1 / rho10 !rho.e(t)
244 mbuf%VAR(i + (k1 + 21 - 1) * nelg) = eint1 * rho1 / rho10 !rho0.e0
245 ELSE
246 mbuf%VAR(i + (k1 + 08 - 1) * nelg) = zero
247 mbuf%VAR(i + (k1 + 21 - 1) * nelg) = zero
248 ENDIF
249 IF (rho20 /= zero) THEN
250 mbuf%VAR(i + (k2 + 08 - 1) * nelg) = eint2 * rho2 / rho20 !rho.e(t)
251 mbuf%VAR(i + (k2 + 21 - 1) * nelg) = eint2 * rho2 / rho20 !rho0.e0
252 ELSE
253 mbuf%VAR(i + (k2 + 08 - 1) * nelg) = zero
254 mbuf%VAR(i + (k2 + 21 - 1) * nelg) = zero
255 ENDIF
256 IF (rho30 /= zero) THEN
257 mbuf%VAR(i + (k3 + 08 - 1) * nelg) = eint3 * rho3 / rho30 !rho.e(t)
258 mbuf%VAR(i + (k3 + 21 - 1) * nelg) = eint3 * rho3 / rho30 !rho0.e0
259 ELSE
260 mbuf%VAR(i + (k3 + 08 - 1) * nelg) = zero
261 mbuf%VAR(i + (k3 + 21 - 1) * nelg) = zero
262 ENDIF
263
264 gbuf%RHO(i) = alpha1 * rho1 + alpha2 * rho2 + alpha3 * rho3 + alpha4 * rho40
265 gbuf%EINT(i) = alpha1 * eint1+ alpha2 * eint2+ alpha3 * eint3+ alpha4 * eint4
266
267 gbuf%SIG(i) = - pgrav
268 gbuf%SIG(i + nelg) = - pgrav
269 gbuf%SIG(i + 2 * nelg) = - pgrav
270
271 mbuf%VAR(i + 3 * nelg) = pgrav
272
273 ENDDO
274
275 END SUBROUTINE inigrav_m51
276
277!||====================================================================
278!|| polysolve ../starter/source/initial_conditions/inigrav/inigrav_m51.F
279!||--- called by ------------------------------------------------------
280!|| inigrav_m51 ../starter/source/initial_conditions/inigrav/inigrav_m51.F
281!||====================================================================
282 SUBROUTINE polysolve(C0, C1, C2, C3, C4, C5, PRES, EINT, MU, RHO0, VOL0, PEXT, UID)
283C-----------------------------------------------
284C I m p l i c i t T y p e s
285C-----------------------------------------------
286#include "implicit_f.inc"
287C-----------------------------------------------
288C D u m m y A r g u m e n t s
289C-----------------------------------------------
290 my_real, INTENT(IN) :: c0, c1, c2, c3, c4, c5, pres, pext
291 my_real, INTENT(OUT) :: mu
292 my_real, INTENT(INOUT) :: rho0, eint, vol0
293 INTEGER, INTENT(IN) :: UID !for debug purpose
294C-----------------------------------------------
295C L o c a l V a r i a b l e s
296C-----------------------------------------------
297 my_real :: tol, error, func, dfunc, incr, vol, vol_prev,pp
298 INTEGER :: ITER, MAX_ITER
299 LOGICAL :: CONT
300C-----------------------------------------------
301C S o u r c e L i n e s
302C-----------------------------------------------
303 iter = 0
304 max_iter = 30
305 tol = em4
306 cont = .true.
307 vol = vol0
308 vol_prev = vol0
309 IF(vol0==zero)vol0=one
310 !Initial guess for MU
311 mu = zero
312 DO WHILE (cont .AND. iter < max_iter)
313 !pressure difference
314 func = (c0 + c1 * mu + c2 * max(zero, mu)*mu + c3 * max(zero, mu)*mu*mu + (c4 + c5 * mu) * eint) - pres
315 !total derivative
316 dfunc = c1 + two * c2 * max(zero, mu) + three * c3 * mu*max(mu,zero) + c5 * eint
317 dfunc = dfunc + (pext+func+pres)/(one+mu)/(one+mu)*(c4+c5*mu)
318 IF(dfunc==zero)EXIT !iform=6 cell with no adj elem (not automatically initialized)
319 ! MU INCREMENT
320 incr = - func / dfunc
321 If(incr==zero)EXIT !P(mu=0)=PRES nothing to do
322 error = abs(func)
323 IF(abs(mu+incr)>em20) error = max(error, abs(incr / (mu + incr)))
324 mu = mu + incr
325 !ENERGY INCREMENT
326 vol = vol0/(one+mu)
327 eint = eint - (pext+func+pres)*(vol-vol_prev)/vol0
328 vol_prev = vol
329 iter = iter + 1
330 IF (error < tol) THEN
331 cont = .false.
332 ENDIF
333 ENDDO
334 !IF(UID==16866)print *,"16866, P=", (C0 + C1 * MU + C2 * MAX(ZERO, MU) ** 2 + C3 * MU ** 3 + (C4 + C5 * MU) * EINT)
335 !IF(UID==16866)print *,"16866, E=", EINT
336 !IF(UID==16866)print *,"16866, MU=", MU
337 !IF(UID==16866)print *, "ITER=", ITER
338 END SUBROUTINE polysolve
my_real
#define my_real
Definition cppsort.cpp:32
alpha2
#define alpha2
Definition eval.h:48
inigrav_m51
subroutine inigrav_m51(nelg, nel, ng, matid, ipm, grav0, depth, pm, bufmat, elbuf_tab, psurf, list, ale_connectivity, nv46, ix, nix, nft, bufmatg, iparg)
Definition inigrav_m51.F:33
polysolve
subroutine polysolve(c0, c1, c2, c3, c4, c5, pres, eint, mu, rho0, vol0, pext, uid)
Definition inigrav_m51.F:283
max
#define max(a, b)
Definition macros.h:21
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
starter
program starter
Definition starter.F:39
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:254