OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
daasolv.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!|| daasolv ../engine/source/fluid/daasolv.F
25!||--- called by ------------------------------------------------------
26!|| flow0 ../engine/source/fluid/flow0.F
27!||--- calls -----------------------------------------------------
28!|| arret ../engine/source/system/arret.F
29!|| finter ../engine/source/tools/curve/finter.F
30!|| finter_smooth ../engine/source/tools/curve/finter_smooth.F
31!||--- uses -----------------------------------------------------
32!|| python_funct_mod ../common_source/modules/python_mod.F90
33!|| sensor_mod ../common_source/modules/sensor_mod.F90
34!||====================================================================
35 SUBROUTINE daasolv(NDIM, NNO, NEL,
36 . IFLOW, IBUF, ELEM, SHELL_GA,
37 . RFLOW, NORMAL, TA, AREAF, COSG, DIST,
38 . MFLE, ACCF, PM, PTI, X, V, A, NPC, TF,
39 . NBGAUGE,LGAUGE, GAUGE, NSENSOR, SENSOR_TAB, IGRV, AGRV,
40 . CBEM ,IPIV , NFUNCT, PYTHON,WFEXT)
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE sensor_mod
45 USE python_funct_mod
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 "param_c.inc"
54!#include "com06_c.inc"
55#include "com08_c.inc"
56C-----------------------------------------------
57C D u m m y A r g u m e n t s
58C-----------------------------------------------
59 INTEGER NDIM, NNO, NEL, NBGAUGE, NSENSOR, NFUNCT
60 INTEGER IFLOW(*), IBUF(*), ELEM(NDIM,*), SHELL_GA(*), NPC(*), LGAUGE(3,*)
61 INTEGER IGRV(NIGRV,*), IPIV(*)
62 my_real X(3,*), V(3,*), A(3,*), TF(*), RFLOW(*)
63 my_real NORMAL(3,*), TA(*), AREAF(*), COSG(*), DIST(*), PM(*), ACCF(*), PTI(*)
64 my_real mfle(nel,*), gauge(llgauge,*),agrv(lfacgrv,*)
65 my_real cbem(nel,*)
66 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) :: SENSOR_TAB
67 TYPE (PYTHON_), INTENT(INOUT) :: PYTHON
68 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
69C-----------------------------------------------
70C L o c a l V a r i a b l e s
71C-----------------------------------------------
72#if defined(MYREAL8) && !defined(WITHOUT_LINALG)
73 INTEGER I, J, K, I1, I2, I3, I4, I5, N1, N2, N3, N4, II, JJ, KK
74 INTEGER IFPRES, JFORM, KFORM, IPRES, IWAVE, INTEGR, FREESURF, AFTERFLOW
75 INTEGER GRAV_ID, IFUNC, ISENS, LENBUF, NNO_L, NEL_L, NELMAX, INFO,ISMOOTH
76 my_real X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X4, Y4, Z4,
77 . X0, Y0, Z0, X12, Y12, Z12, X13, Y13, Z13, X24, Y24, Z24,
78 . nrx, nry, nrz, area2, wf(2), wi(4,2),
79 . pm1, dsc, vx0, vy0, vz0
80 my_real rho, ssp, rhoc, tta, dydx, ff, sfpres, dppdt, pmax, theta, norm, fac1, ts, gravity
81 my_real dirwix, dirwiy, dirwiz, dirwrx, dirwry, dirwrz
82 my_real off(nel), pp(nel), vi(nel), ps(nel), fel(nel), pr(nel), vr(nel)
83 my_real dpidt(nel), dprdt(nel), dpmdt(nel), hh(nel)
84 my_real xa, ya, za, fcx, fcy, pmin, ptot
85 my_real apmax, atheta, ratio
86
87 my_real finter,finter_smooth
88 EXTERNAL finter,finter_smooth
89 my_real vect(nel), vect1(nel), vel(nel)
90 my_real xl(3,nno), vl(3,nno)
91
92C Coordonnees et vitesses locales
93 DO i=1,nno
94 ii=ibuf(i)
95 xl(1,i)=x(1,ii)
96 xl(2,i)=x(2,ii)
97 xl(3,i)=x(3,ii)
98 vl(1,i)=v(1,ii)
99 vl(2,i)=v(2,ii)
100 vl(3,i)=v(3,ii)
101 ENDDO
102
103 jform = iflow(4)
104 ipres = iflow(21)
105 iwave = iflow(22)
106 kform = iflow(23)
107 integr = iflow(24)
108 freesurf = iflow(25)
109 afterflow = iflow(26)
110 grav_id = iflow(27)
111 nelmax = iflow(28)
112 rhoc = rflow(1)
113 ssp = rflow(2)
114 rho = rflow(5)
115 fac1 = rflow(8)
116 dsc = rflow(12)
117 xa = rflow(16)
118 ya = rflow(17)
119 za = rflow(18)
120 pmin = rflow(22)
121 apmax = rflow(23)
122 atheta = rflow(24)
123
124 IF(jform == 1) THEN
125 DO i = 1,nel
126 i1 = elem(1,i)
127 i2 = elem(2,i)
128 i3 = elem(3,i)
129 x1 = xl(1,i1)
130 x2 = xl(1,i2)
131 x3 = xl(1,i3)
132 y1 = xl(2,i1)
133 y2 = xl(2,i2)
134 y3 = xl(2,i3)
135 z1 = xl(3,i1)
136 z2 = xl(3,i2)
137 z3 = xl(3,i3)
138 x12= x2-x1
139 y12= y2-y1
140 z12= z2-z1
141 x13= x3-x1
142 y13= y3-y1
143 z13= z3-z1
144 nrx= y12*z13-z12*y13
145 nry= z12*x13-x12*z13
146 nrz= x12*y13-y12*x13
147 area2 = sqrt(nrx**2+nry**2+nrz**2)
148 normal(1,i) = nrx/area2
149 normal(2,i) = nry/area2
150 normal(3,i) = nrz/area2
151 areaf(i) = half*area2
152 IF(iwave == 1)THEN
153C Spherical wave
154 x0=third*x1+third*x2+third*x3
155 y0=third*y1+third*y2+third*y3
156 z0=third*z1+third*z2+third*z3
157 dirwix=x0-rflow(9)
158 dirwiy=y0-rflow(10)
159 dirwiz=z0-rflow(11)
160 norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
161 dirwix=dirwix/norm
162 dirwiy=dirwiy/norm
163 dirwiz=dirwiz/norm
164 IF(freesurf == 2) THEN
165 dirwrx=x0-rflow(13)
166 dirwry=y0-rflow(14)
167 dirwrz=z0-rflow(15)
168 norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
169 dirwrx=dirwrx/norm
170 dirwry=dirwry/norm
171 dirwrz=dirwrz/norm
172 cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
173 ENDIF
174C hydrostatic pressure
175 hh(i)=zero
176 IF(grav_id > 0) THEN
177 hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
178 ENDIF
179 ELSEIF(iwave == 2) THEN
180C Onde plane
181 dirwix=rflow(9)
182 dirwiy=rflow(10)
183 dirwiz=rflow(11)
184C hydrostatic pressure
185 hh(i)=zero
186 ENDIF
187 cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
188C vitesse N-1/2
189 vx0 = third*vl(1,i1)+third*vl(1,i2)+third*vl(1,i3)
190 vy0 = third*vl(2,i1)+third*vl(2,i2)+third*vl(2,i3)
191 vz0 = third*vl(3,i1)+third*vl(3,i2)+third*vl(3,i3)
192 vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
193 ENDDO
194 ELSEIF(jform == 2) THEN
195 wi(1,1)=fourth
196 wi(2,1)=fourth
197 wi(3,1)=fourth
198 wi(4,1)=fourth
199 wi(1,2)=third
200 wi(2,2)=third
201 wi(3,2)=one_over_6
202 wi(4,2)=one_over_6
203 DO i = 1,nel
204 i1 = elem(1,i)
205 i2 = elem(2,i)
206 i3 = elem(3,i)
207 i4 = elem(4,i)
208 i5 = elem(5,i)
209 x1 = xl(1,i1)
210 x2 = xl(1,i2)
211 x3 = xl(1,i3)
212 x4 = xl(1,i4)
213 y1 = xl(2,i1)
214 y2 = xl(2,i2)
215 y3 = xl(2,i3)
216 y4 = xl(2,i4)
217 z1 = xl(3,i1)
218 z2 = xl(3,i2)
219 z3 = xl(3,i3)
220 z4 = xl(3,i4)
221 x13= x3-x1
222 y13= y3-y1
223 z13= z3-z1
224 x24= x4-x2
225 y24= y4-y2
226 z24= z4-z2
227 nrx=y13*z24-z13*y24
228 nry=z13*x24-x13*z24
229 nrz=x13*y24-y13*x24
230 area2 = sqrt(nrx**2+nry**2+nrz**2)
231 normal(1,i) = nrx/area2
232 normal(2,i) = nry/area2
233 normal(3,i) = nrz/area2
234 areaf(i) = half*area2
235 IF(iwave == 1)THEN
236C Spherical wave
237 x0=wi(1,i5)*x1+wi(2,i5)*x2+wi(3,i5)*x3+wi(4,i5)*x4
238 y0=wi(1,i5)*y1+wi(2,i5)*y2+wi(3,i5)*y3+wi(4,i5)*y4
239 z0=wi(1,i5)*z1+wi(2,i5)*z2+wi(3,i5)*z3+wi(4,i5)*z4
240 dirwix=x0-rflow(9)
241 dirwiy=y0-rflow(10)
242 dirwiz=z0-rflow(11)
243 norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
244 dirwix=dirwix/norm
245 dirwiy=dirwiy/norm
246 dirwiz=dirwiz/norm
247 IF(freesurf == 2) THEN
248 dirwrx=x0-rflow(13)
249 dirwry=y0-rflow(14)
250 dirwrz=z0-rflow(15)
251 norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
252 dirwrx=dirwrx/norm
253 dirwry=dirwry/norm
254 dirwrz=dirwrz/norm
255 cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
256 ENDIF
257C hydrostatic pressure
258 hh(i)=zero
259 IF(grav_id > 0) THEN
260 hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
261 ENDIF
262 ELSEIF(iwave == 2) THEN
263C Onde plane
264 dirwix=rflow(9)
265 dirwiy=rflow(10)
266 dirwiz=rflow(11)
267C hydrostatic pressure
268 hh(i)=zero
269 ENDIF
270 cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
271C vitesse N-1/2
272 vx0 = wi(1,i5)*vl(1,i1)+wi(2,i5)*vl(1,i2)+wi(3,i5)*vl(1,i3)+wi(4,i5)*vl(1,i4)
273 vy0 = wi(1,i5)*vl(2,i1)+wi(2,i5)*vl(2,i2)+wi(3,i5)*vl(2,i3)+wi(4,i5)*vl(2,i4)
274 vz0 = wi(1,i5)*vl(3,i1)+wi(2,i5)*vl(3,i2)+wi(3,i5)*vl(3,i3)+wi(4,i5)*vl(3,i4)
275 vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
276 ENDDO
277
278 ENDIF
279C-----------------------------------------------
280C 1. Pression incidente
281C-----------------------------------------------
282 DO i=1,nel
283 off(i) = zero
284 tta = tt-ta(i)
285 IF(tta > zero) off(i) = one
286 ENDDO
287 IF(ipres == 1) THEN
288 IF(iwave == 1)THEN
289 DO i=1,nel
290 tta = tt-ta(i)
291 ratio = dsc/dist(i)
292 pmax = rflow(6)*ratio**apmax
293 theta = rflow(7)*ratio**atheta
294 pp(i) = pmax*exp(-tta/theta)*off(i)
295 dpidt(i) = -pp(i)*cosg(i)/theta
296 ENDDO
297 ELSEIF(iwave == 2) THEN
298 pmax = rflow(6)
299 theta = rflow(7)
300 DO i=1,nel
301 tta = tt-ta(i)
302 pp(i) = pmax*exp(-tta/theta)*off(i)
303 dpidt(i) = -pp(i)*cosg(i)/theta
304 ENDDO
305 ENDIF
306 ELSEIF(ipres == 2) THEN
307 ifpres = iflow(7)
308 sfpres = rflow(3)
309 IF(iwave == 1)THEN
310 DO i=1,nel
311 tta = tt-ta(i)
312 ratio = dsc/dist(i)
313 IF(ifpres > 0) THEN
314 ismooth = npc(2*nfunct+ifpres+1)
315!! PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
316 IF (ismooth == 0) THEN
317 pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
318 ELSE IF(ismooth > 0) THEN
319 pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
320 ELSE
321 ismooth = -ismooth ! function id
322 CALL python_call_funct1d(python, ismooth,tta,pp(i))
323 pp(i) = pp(i)*sfpres*ratio*off(i)
324 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
325 ENDIF ! IF (ISMOOTH == 0)
326 dpidt(i) = dppdt*cosg(i)*ratio*off(i)
327 ELSE
328 pp(i) = sfpres*ratio*off(i)
329 dpidt(i) = zero
330 ENDIF
331 ENDDO
332 ELSEIF(iwave == 2) THEN
333 DO i=1,nel
334 tta = tt-ta(i)
335 IF(ifpres > 0) THEN
336 ismooth = npc(2*nfunct+ifpres+1)
337!! PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*OFF(I)
338 IF (ismooth == 0) THEN
339 pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*off(i)
340 ELSE IF (ismooth > 0) THEN
341 pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*off(i)
342 ELSE
343 ismooth = -ismooth ! function id
344 CALL python_call_funct1d(python, ismooth,tta,pp(i))
345 pp(i) = pp(i)*sfpres*off(i)
346 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
347 ENDIF ! IF (ISMOOTH == 0)
348 dpidt(i) = dppdt*cosg(i)*off(i)
349 ELSE
350 pp(i) = sfpres*off(i)
351 dpidt(i) = zero
352 ENDIF
353 ENDDO
354 ENDIF
355 ENDIF
356C-----------------------------------------------
357C Vitesse onde incidente
358C-----------------------------------------------
359 DO i=1,nel
360 vi(i)=pp(i)*cosg(i)/rhoc
361 ENDDO
362 IF(afterflow == 2) THEN
363 IF(iwave == 1)THEN
364C add afterflow velocity, no afterflow velocity for plane wave
365 IF(ipres == 1) THEN
366 DO i=1,nel
367 ratio = dsc/dist(i)
368 pmax = rflow(6)*ratio**apmax
369 theta = rflow(7)*ratio**atheta
370 vi(i) = vi(i) + (pmax-pp(i))*cosg(i)*theta/(rho*dist(i))
371 ENDDO
372 ELSEIF(ipres == 2) THEN
373 DO i=1,nel
374 pti(i) = pti(i) + pp(i)*dt1
375 vi(i) = vi(i) + pti(i)*cosg(i)/(rho*dist(i))
376 ENDDO
377 ENDIF
378 ENDIF
379 ENDIF
380C-----------------------------------------------
381C 2. Pression Onde Reflechie
382C-----------------------------------------------
383 IF(freesurf == 2) THEN
384 DO i=1,nel
385 off(i) = zero
386 tta = tt-ta(i+nel)
387 IF(tta > zero) off(i) = one
388 ENDDO
389 IF(ipres == 1) THEN
390 DO i=1,nel
391 j = i+nel
392 tta = tt-ta(j)
393 ratio = dsc/dist(j)
394 pmax = rflow(6)*ratio**apmax
395 theta = rflow(7)*ratio**atheta
396 pr(i) = -pmax*exp(-tta/theta)*off(i)
397 dprdt(i) = -pr(i)*cosg(j)/theta
398 ENDDO
399 ELSEIF(ipres == 2) THEN
400 DO i=1,nel
401 j = i+nel
402 tta = tt-ta(j)
403 ratio = dsc/dist(j)
404 IF(ifpres > 0) THEN
405 ismooth = npc(2*nfunct+ifpres+1)
406!! PR(I) = -SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
407 IF (ismooth == 0) THEN
408 pr(i) = -sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
409 ELSE IF (ismooth > 0) THEN
410 pr(i) = -sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
411 ELSE
412 ismooth = -ismooth ! function id
413 CALL python_call_funct1d(python, ismooth,tta,pr(i))
414 pr(i) = -pr(i)*sfpres*ratio*off(i)
415 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
416 ENDIF ! IF (ISMOOTH == 0)
417 dprdt(i) = dppdt*cosg(j)*off(i)
418 ELSE
419 pr(i) = -sfpres*ratio*off(i)
420 dprdt(i) = zero
421 ENDIF
422 ENDDO
423 ENDIF
424C-----------------------------------------------
425C Vitesse onde reflechie
426C-----------------------------------------------
427 DO i=1,nel
428 vr(i)=pr(i)*cosg(i+nel)/rhoc
429 ENDDO
430 ELSE ! pas de surface libre
431 DO i=1,nel
432 pr(i) = zero
433 vr(i) = zero
434 dprdt(i) = zero
435 ENDDO
436 ENDIF
437C-----------------------------------------------
438C 3. Scattered pressure Schema Euler ordre 1
439C-----------------------------------------------
440 IF(kform == 1) THEN
441C DAA-1
442 IF(nel < nelmax) THEN
443 DO i=1,nel
444 dpmdt(i) = rhoc*accf(i)
445 vect(i) = rhoc*areaf(i)*(pm(i)-rhoc*(vi(i)+vr(i)))
446 ENDDO
447 vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
448 DO i=1,nel
449 dpmdt(i) = dpmdt(i) - vect(i)
450 ENDDO
451 IF(integr == 1) THEN
452 DO i=1,nel
453 pm(i) = pm(i) + dpmdt(i)*dt1
454 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
455 ENDDO
456
457 ELSEIF(integr == 2) THEN
458C Prediction
459 DO i=1,nel
460 ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
461 ENDDO
462C Correction
463 DO i=1,nel
464 dpmdt(i) = rhoc*accf(i)
465 vect(i) = rhoc*areaf(i)*ps(i)
466 ENDDO
467 vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
468 DO i=1,nel
469 dpmdt(i) = dpmdt(i) - vect(i)
470 ENDDO
471 DO i=1,nel
472 pm(i) = pm(i) + dpmdt(i)*dt1
473 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
474 ENDDO
475 ENDIF
476C
477 ELSE ! resolution systeme MFLE*X=Y
478 IF(dt1 == zero) THEN
479C L U decomposition
480 CALL dgetrf(nel, nel, mfle, nel, ipiv, info)
481 ENDIF
482 DO i=1,nel
483 vect(i) = two*ssp*(rhoc*(vi(i)+vr(i))-pm(i))
484 ENDDO
485 CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
486C Resolution
487 CALL dgetrs('N', nel, 1, mfle, nel, ipiv, vect1, nel, info)
488 vect(1:nel) = matmul(cbem(1:nel,1:nel), vect1(1:nel))
489 DO i=1,nel
490 IF(pp(i) /= zero) THEN
491 dpmdt(i) = rhoc*accf(i) + vect(i)
492 ELSE
493 dpmdt(i)=zero
494 ENDIF
495 ENDDO
496 IF(integr == 1) THEN
497 DO i=1,nel
498 pm(i) = pm(i) + dpmdt(i)*dt1
499 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
500 ENDDO
501
502 ELSEIF(integr == 2) THEN
503C Prediction
504 DO i=1,nel
505 ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
506 ENDDO
507C Correction
508 DO i=1,nel
509 vect(i) = -two*ssp*ps(i)
510 ENDDO
511 CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
512C Resolution
513 CALL dgetrs('N', nel, 1, mfle, nel, ipiv, vect1, nel, info)
514 vect(1:nel) = matmul(cbem(1:nel,1:nel), vect1(1:nel))
515 DO i=1,nel
516 dpmdt(i) = rhoc*accf(i) + vect(i)
517 ENDDO
518 DO i=1,nel
519 pm(i) = pm(i) + dpmdt(i)*dt1
520 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
521 ENDDO
522 ENDIF
523 ENDIF
524C
525 ELSEIF(kform == 2) THEN
526C High Frequency Approximation
527 DO i=1,nel
528 dpmdt(i) = rhoc*accf(i)
529 ENDDO
530 DO i=1,nel
531 pm(i) = pm(i) + dpmdt(i)*dt1
532 ps(i) = pm(i) - rhoc*vi(i)
533 ENDDO
534C
535 ELSEIF(kform == 3) THEN
536C Virtual Mass Approximation
537 DO i=1,nel
538 ps(i) = zero
539 DO j=1,nel
540 ps(i) = ps(i) + mfle(i,j)*(accf(j)-(dpidt(j)+dprdt(j))/rhoc)
541 ENDDO
542 ps(i) = ps(i)/areaf(i)
543 ENDDO
544 ENDIF
545C-----------------------------------------------
546C 4. gravity
547C-----------------------------------------------
548 gravity = zero
549 IF(grav_id > 0) THEN
550 fcy = agrv(1,grav_id)
551 fcx = agrv(2,grav_id)
552 ifunc = igrv(3,grav_id)
553 isens = 0
554 DO k=1,nsensor
555 IF(igrv(6,grav_id) == sensor_tab(k)%SENS_ID) isens=k ! do it in starter !!!
556 ENDDO
557 IF(isens==0)THEN
558 ts = tt
559 ELSE
560 ts = tt-sensor_tab(isens)%TSTART
561 ENDIF
562 ismooth = 0
563 IF (ifunc > 0) THEN
564 ismooth = npc(2*nfunct+ifunc+1)
565!! GRAVITY = FCY*FINTER(IFUNC,TS*FCX,NPC,TF,DYDX)
566 IF (ismooth == 0) THEN
567 gravity = fcy*finter(ifunc,ts*fcx,npc,tf,dydx)
568 ELSE IF(ismooth > 0) THEN
569 gravity = fcy*finter_smooth(ifunc,ts*fcx,npc,tf,dydx)
570 ELSE
571 ismooth = -ismooth ! function id
572 CALL python_call_funct1d(python, ismooth,ts*fcx,gravity)
573 gravity = gravity*fcy
574 ENDIF ! IF (ISMOOTH = 0)
575 ELSE
576 gravity = fcy
577 ENDIF
578 ENDIF
579C-----------------------------------------------
580C 5. compute total pressure forces
581C-----------------------------------------------
582 DO i=1,nel
583 ptot = pp(i)+ps(i)+pr(i)+rho*hh(i)*abs(gravity)
584 fel(i) = areaf(i)*max(ptot,pmin)
585 ENDDO
586C
587 IF(jform == 1) THEN
588 DO i=1,nel
589 i1 = elem(1,i)
590 i2 = elem(2,i)
591 i3 = elem(3,i)
592 n1 = ibuf(i1)
593 n2 = ibuf(i2)
594 n3 = ibuf(i3)
595 nrx= normal(1,i)
596 nry= normal(2,i)
597 nrz= normal(3,i)
598 ff = fac1*third*fel(i)
599 wfext = wfext-fel(i)*vel(i)*dt1
600 a(1,n1) = a(1,n1) + ff*nrx
601 a(2,n1) = a(2,n1) + ff*nry
602 a(3,n1) = a(3,n1) + ff*nrz
603 a(1,n2) = a(1,n2) + ff*nrx
604 a(2,n2) = a(2,n2) + ff*nry
605 a(3,n2) = a(3,n2) + ff*nrz
606 a(1,n3) = a(1,n3) + ff*nrx
607 a(2,n3) = a(2,n3) + ff*nry
608 a(3,n3) = a(3,n3) + ff*nrz
609 ENDDO
610 ELSEIF(jform == 2) THEN
611 wf(1)=fac1*fourth
612 wf(2)=fac1*third
613 DO i=1,nel
614 i1 = elem(1,i)
615 i2 = elem(2,i)
616 i3 = elem(3,i)
617 i4 = elem(4,i)
618 i5 = elem(5,i)
619 n1 = ibuf(i1)
620 n2 = ibuf(i2)
621 n3 = ibuf(i3)
622 n4 = ibuf(i4)
623 nrx= normal(1,i)
624 nry= normal(2,i)
625 nrz= normal(3,i)
626 ff = wf(i5)*fel(i)
627 wfext = wfext-fel(i)*vel(i)*dt1
628 a(1,n1) = a(1,n1) + ff*nrx
629 a(2,n1) = a(2,n1) + ff*nry
630 a(3,n1) = a(3,n1) + ff*nrz
631 a(1,n2) = a(1,n2) + ff*nrx
632 a(2,n2) = a(2,n2) + ff*nry
633 a(3,n2) = a(3,n2) + ff*nrz
634 a(1,n3) = a(1,n3) + ff*nrx
635 a(2,n3) = a(2,n3) + ff*nry
636 a(3,n3) = a(3,n3) + ff*nrz
637 IF(i5 == 2) cycle
638 a(1,n4) = a(1,n4) + ff*nrx
639 a(2,n4) = a(2,n4) + ff*nry
640 a(3,n4) = a(3,n4) + ff*nrz
641 ENDDO
642 ENDIF
643C-----------------------------------------------------------
644C 5. Gauge
645C-----------------------------------------------------------
646 IF(jform == 2) THEN
647 DO i=1,nbgauge
648 gauge(30,i)=zero
649 IF(lgauge(1,i) == 0 .AND. lgauge(3,i) < 0 ) THEN
650 i1=shell_ga(i)
651 IF(i1 > 0) THEN
652 ptot=pp(i1)+ps(i1)+pr(i1)+rho*hh(i1)*abs(gravity)
653 gauge(30,i)=max(ptot,pmin)
654 ENDIF
655 ENDIF
656 ENDDO
657 ENDIF
658#else
659 CALL arret(5)
660#endif
661
662 RETURN
663 END
664
daasolv
subroutine daasolv(ndim, nno, nel, iflow, ibuf, elem, shell_ga, rflow, normal, ta, areaf, cosg, dist, mfle, accf, pm, pti, x, v, a, npc, tf, nbgauge, lgauge, gauge, nsensor, sensor_tab, igrv, agrv, cbem, ipiv, nfunct, python, wfext)
Definition daasolv.F:41
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
dgetrs
subroutine dgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
DGETRS
Definition dgetrs.f:121
dgetrf
subroutine dgetrf(m, n, a, lda, ipiv, info)
DGETRF
Definition dgetrf.f:108
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
max
#define max(a, b)
Definition macros.h:21
arret
subroutine arret(nn)
Definition arret.F:87