OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
daasolvp.F File Reference
#include "implicit_f.inc"
#include "param_c.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine daasolvp (ndim, nno, nel, nel_loc, iflow, ibuf, elem, ibufl, shell_ga, cnp, 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_l, mfle_l, ibufelr, ibufelc, nfunct, python, wfext)

Function/Subroutine Documentation

◆ daasolvp()

subroutine daasolvp ( integer ndim,
integer nno,
integer nel,
integer nel_loc,
integer, dimension(*) iflow,
integer, dimension(*) ibuf,
integer, dimension(ndim,*) elem,
integer, dimension(*) ibufl,
integer, dimension(*) shell_ga,
integer, dimension(*) cnp,
rflow,
normal,
ta,
areaf,
cosg,
dist,
mfle,
accf,
pm,
pti,
x,
v,
a,
integer, dimension(*) npc,
tf,
integer nbgauge,
integer, dimension(3,*) lgauge,
gauge,
integer nsensor,
type (sensor_str_), dimension(nsensor) sensor_tab,
integer, dimension(nigrv,*) igrv,
agrv,
cbem,
integer, dimension(*) ipiv_l,
mfle_l,
integer, dimension(nel) ibufelr,
integer, dimension(nel) ibufelc,
integer nfunct,
type (python_), intent(inout) python,
double precision, intent(inout) wfext )

Definition at line 35 of file daasolvp.F.

42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE sensor_mod
46 USE python_funct_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "param_c.inc"
55!#include "com06_c.inc"
56#include "com08_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER NDIM, NNO, NEL, NBGAUGE, NSENSOR, NEL_LOC,NFUNCT
61 INTEGER IFLOW(*), IBUF(*), ELEM(NDIM,*), SHELL_GA(*), NPC(*), LGAUGE(3,*)
62 INTEGER IBUFL(*), CNP(*), IGRV(NIGRV,*)
63 INTEGER IPIV_L(*), IBUFELR(NEL), IBUFELC(NEL)
64 my_real x(3,*), v(3,*), a(3,*), tf(*), rflow(*)
65 my_real normal(3,*), ta(*), areaf(*), cosg(*), dist(*), pm(*), accf(*), pti(*)
66 my_real mfle(nel,*), gauge(llgauge,*), agrv(lfacgrv,*)
67 my_real cbem(nel,*), mfle_l(nel_loc,*)
68 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) :: SENSOR_TAB
69 TYPE (PYTHON_), INTENT(INOUT) :: PYTHON
70 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER I, J, K, L, I1, I2, I3, I4, I5, N1, N2, N3, N4, II, JJ, KK
75 INTEGER IFPRES, JFORM, KFORM, IPRES, IWAVE, INTEGR, FREESURF, AFTERFLOW
76 INTEGER GRAV_ID, IFUNC, ISENS, LENBUF, NNO_L, NEL_L, NELMAX, INFO,ISMOOTH
77 INTEGER NEBUF(NEL)
78 INTEGER NBLOC, NPROW, NPCOL, ICTXT, DESCA(9), DESCB(9), MYROW, MYCOL, MXLLDM, NUMROC
79 EXTERNAL numroc
80 my_real x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4,
81 . x0, y0, z0, x12, y12, z12, x13, y13, z13, x24, y24, z24,
82 . nrx, nry, nrz, area2, wf(2), wi(4,2),
83 . pm1, dsc, vx0, vy0, vz0
84 my_real rho, ssp, rhoc, tta, dydx, ff, sfpres, dppdt, pmax, theta, norm, fac1, ts, gravity
85 my_real dirwix, dirwiy, dirwiz, dirwrx, dirwry, dirwrz
86 my_real off(nel), pp(nel), vi(nel), ps(nel), fel(nel), pr(nel), vr(nel)
87 my_real dpidt(nel), dprdt(nel), dpmdt(nel), hh(nel)
88 my_real xa, ya, za, fcx, fcy, pmin, ptot
89 my_real apmax, atheta, ratio
90 my_real finter,finter_smooth
91 EXTERNAL finter,finter_smooth
92 my_real vect(nel), vect1(nel), vel(nel)
93 my_real xl(3,nno), vl(3,nno)
94 my_real, ALLOCATABLE :: sbuf(:), rbuf(:)
95
96#if defined(MPI) && defined(MYREAL8) && !defined(WITHOUT_LINALG)
97C Coordonnees et vitesses locales SPMD
98 nno_l=iflow(16)
99 lenbuf=6*nno
100 ALLOCATE(sbuf(lenbuf), rbuf(lenbuf))
101 sbuf(1:lenbuf)=zero
102 rbuf(1:lenbuf)=zero
103 DO i=1,nno_l
104 ii=ibufl(i)
105 jj=ibuf(ii)
106 kk=6*(ii-1)
107 sbuf(kk+1)=x(1,jj)/cnp(ii)
108 sbuf(kk+2)=x(2,jj)/cnp(ii)
109 sbuf(kk+3)=x(3,jj)/cnp(ii)
110 sbuf(kk+4)=v(1,jj)/cnp(ii)
111 sbuf(kk+5)=v(2,jj)/cnp(ii)
112 sbuf(kk+6)=v(3,jj)/cnp(ii)
113 ENDDO
114
115 CALL spmd_fl_sum(sbuf, lenbuf, rbuf)
116
117 DO i=1,nno
118 k=6*(i-1)
119 xl(1,i)=rbuf(k+1)
120 xl(2,i)=rbuf(k+2)
121 xl(3,i)=rbuf(k+3)
122 vl(1,i)=rbuf(k+4)
123 vl(2,i)=rbuf(k+5)
124 vl(3,i)=rbuf(k+6)
125 ENDDO
126 DEALLOCATE(sbuf, rbuf)
127
128 jform = iflow(4)
129 nbloc = iflow(12)
130 nprow = iflow(18)
131 npcol = iflow(19)
132 ipres = iflow(21)
133 iwave = iflow(22)
134 kform = iflow(23)
135 integr = iflow(24)
136 freesurf = iflow(25)
137 afterflow = iflow(26)
138 grav_id = iflow(27)
139 nelmax = iflow(28)
140 rhoc = rflow(1)
141 ssp = rflow(2)
142 rho = rflow(5)
143 fac1 = rflow(8)
144 dsc = rflow(12)
145 xa = rflow(16)
146 ya = rflow(17)
147 za = rflow(18)
148 pmin = rflow(22)
149 apmax = rflow(23)
150 atheta = rflow(24)
151
152 IF(jform == 1) THEN
153 DO i = 1,nel
154 i1 = elem(1,i)
155 i2 = elem(2,i)
156 i3 = elem(3,i)
157 x1 = xl(1,i1)
158 x2 = xl(1,i2)
159 x3 = xl(1,i3)
160 y1 = xl(2,i1)
161 y2 = xl(2,i2)
162 y3 = xl(2,i3)
163 z1 = xl(3,i1)
164 z2 = xl(3,i2)
165 z3 = xl(3,i3)
166 x12= x2-x1
167 y12= y2-y1
168 z12= z2-z1
169 x13= x3-x1
170 y13= y3-y1
171 z13= z3-z1
172 nrx= y12*z13-z12*y13
173 nry= z12*x13-x12*z13
174 nrz= x12*y13-y12*x13
175 area2 = sqrt(nrx**2+nry**2+nrz**2)
176 normal(1,i) = nrx/area2
177 normal(2,i) = nry/area2
178 normal(3,i) = nrz/area2
179 areaf(i) = half*area2
180 IF(iwave == 1)THEN
181C Spherical wave
182 x0=third*x1+third*x2+third*x3
183 y0=third*y1+third*y2+third*y3
184 z0=third*z1+third*z2+third*z3
185 dirwix=x0-rflow(9)
186 dirwiy=y0-rflow(10)
187 dirwiz=z0-rflow(11)
188 norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
189 dirwix=dirwix/norm
190 dirwiy=dirwiy/norm
191 dirwiz=dirwiz/norm
192 IF(freesurf == 2) THEN
193 dirwrx=x0-rflow(13)
194 dirwry=y0-rflow(14)
195 dirwrz=z0-rflow(15)
196 norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
197 dirwrx=dirwrx/norm
198 dirwry=dirwry/norm
199 dirwrz=dirwrz/norm
200 cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
201 ENDIF
202C hydrostatic pressure
203 hh(i)=zero
204 IF(grav_id > 0) THEN
205 hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
206 ENDIF
207 ELSEIF(iwave == 2) THEN
208C Onde plane
209 dirwix=rflow(9)
210 dirwiy=rflow(10)
211 dirwiz=rflow(11)
212C hydrostatic pressure
213 hh(i)=zero
214 ENDIF
215 cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
216C vitesse N-1/2
217 vx0 = third*vl(1,i1)+third*vl(1,i2)+third*vl(1,i3)
218 vy0 = third*vl(2,i1)+third*vl(2,i2)+third*vl(2,i3)
219 vz0 = third*vl(3,i1)+third*vl(3,i2)+third*vl(3,i3)
220 vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
221 ENDDO
222 ELSEIF(jform == 2) THEN
223 wi(1,1)=fourth
224 wi(2,1)=fourth
225 wi(3,1)=fourth
226 wi(4,1)=fourth
227 wi(1,2)=third
228 wi(2,2)=third
229 wi(3,2)=one_over_6
230 wi(4,2)=one_over_6
231 DO i = 1,nel
232 i1 = elem(1,i)
233 i2 = elem(2,i)
234 i3 = elem(3,i)
235 i4 = elem(4,i)
236 i5 = elem(5,i)
237 x1 = xl(1,i1)
238 x2 = xl(1,i2)
239 x3 = xl(1,i3)
240 x4 = xl(1,i4)
241 y1 = xl(2,i1)
242 y2 = xl(2,i2)
243 y3 = xl(2,i3)
244 y4 = xl(2,i4)
245 z1 = xl(3,i1)
246 z2 = xl(3,i2)
247 z3 = xl(3,i3)
248 z4 = xl(3,i4)
249 x13= x3-x1
250 y13= y3-y1
251 z13= z3-z1
252 x24= x4-x2
253 y24= y4-y2
254 z24= z4-z2
255 nrx=y13*z24-z13*y24
256 nry=z13*x24-x13*z24
257 nrz=x13*y24-y13*x24
258 area2 = sqrt(nrx**2+nry**2+nrz**2)
259 normal(1,i) = nrx/area2
260 normal(2,i) = nry/area2
261 normal(3,i) = nrz/area2
262 areaf(i) = half*area2
263 IF(iwave == 1)THEN
264C Spherical wave
265 x0=wi(1,i5)*x1+wi(2,i5)*x2+wi(3,i5)*x3+wi(4,i5)*x4
266 y0=wi(1,i5)*y1+wi(2,i5)*y2+wi(3,i5)*y3+wi(4,i5)*y4
267 z0=wi(1,i5)*z1+wi(2,i5)*z2+wi(3,i5)*z3+wi(4,i5)*z4
268 dirwix=x0-rflow(9)
269 dirwiy=y0-rflow(10)
270 dirwiz=z0-rflow(11)
271 norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
272 dirwix=dirwix/norm
273 dirwiy=dirwiy/norm
274 dirwiz=dirwiz/norm
275 IF(freesurf == 2) THEN
276 dirwrx=x0-rflow(13)
277 dirwry=y0-rflow(14)
278 dirwrz=z0-rflow(15)
279 norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
280 dirwrx=dirwrx/norm
281 dirwry=dirwry/norm
282 dirwrz=dirwrz/norm
283 cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
284 ENDIF
285C hydrostatic pressure
286 hh(i)=zero
287 IF(grav_id > 0) THEN
288 hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
289 ENDIF
290 ELSEIF(iwave == 2) THEN
291C Onde plane
292 dirwix=rflow(9)
293 dirwiy=rflow(10)
294 dirwiz=rflow(11)
295C hydrostatic pressure
296 hh(i)=zero
297 ENDIF
298 cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
299C vitesse N-1/2
300 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)
301 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)
302 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)
303 vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
304 ENDDO
305
306 ENDIF
307C-----------------------------------------------
308C 1. Pression incidente
309C-----------------------------------------------
310 DO i=1,nel
311 off(i) = zero
312 tta = tt-ta(i)
313 IF(tta > zero) off(i) = one
314 ENDDO
315 IF(ipres == 1) THEN
316 IF(iwave == 1)THEN
317 DO i=1,nel
318 tta = tt-ta(i)
319 ratio = dsc/dist(i)
320 pmax = rflow(6)*ratio**apmax
321 theta = rflow(7)*ratio**atheta
322 pp(i) = pmax*exp(-tta/theta)*off(i)
323 dpidt(i) = -pp(i)*cosg(i)/theta
324 ENDDO
325 ELSEIF(iwave == 2) THEN
326 pmax = rflow(6)
327 theta = rflow(7)
328 DO i=1,nel
329 tta = tt-ta(i)
330 pp(i) = pmax*exp(-tta/theta)*off(i)
331 dpidt(i) = -pp(i)*cosg(i)/theta
332 ENDDO
333 ENDIF
334 ELSEIF(ipres == 2) THEN
335 ifpres = iflow(7)
336 sfpres = rflow(3)
337 IF(iwave == 1)THEN
338 DO i=1,nel
339 tta = tt-ta(i)
340 ratio = dsc/dist(i)
341 IF(ifpres > 0) THEN
342 ismooth = npc(2*nfunct+ifpres+1)
343!! PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
344 IF (ismooth == 0) THEN
345 pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
346 ELSE IF(ismooth > 0) THEN
347 pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
348 ELSE
349 ismooth = -ismooth ! function id
350 CALL python_call_funct1d(python, ismooth,tta,pp(i))
351 pp(i) = pp(i)*sfpres*ratio*off(i)
352 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
353 ENDIF ! IF (ISMOOTH == 0)
354 dpidt(i) = dppdt*cosg(i)*ratio*off(i)
355 ELSE
356 pp(i) = sfpres*ratio*off(i)
357 dpidt(i) = zero
358 ENDIF
359 ENDDO
360 ELSEIF(iwave == 2) THEN
361 DO i=1,nel
362 tta = tt-ta(i)
363 IF(ifpres > 0) THEN
364 ismooth = npc(2*nfunct+ifpres+1)
365!! PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*OFF(I)
366 IF (ismooth == 0) THEN
367 pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*off(i)
368 ELSE IF (ismooth > 0) THEN
369 pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*off(i)
370 ELSE
371 ismooth = -ismooth ! function id
372 CALL python_call_funct1d(python, ismooth,tta,pp(i))
373 pp(i) = pp(i)*sfpres*off(i)
374 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
375 ENDIF ! IF (ISMOOTH == 0)
376 dpidt(i) = dppdt*cosg(i)*off(i)
377 ELSE
378 pp(i) = sfpres*off(i)
379 dpidt(i) = zero
380 ENDIF
381 ENDDO
382 ENDIF
383 ENDIF
384C-----------------------------------------------
385C Vitesse onde incidente
386C-----------------------------------------------
387 DO i=1,nel
388 vi(i)=pp(i)*cosg(i)/rhoc
389 ENDDO
390 IF(afterflow == 2) THEN
391 IF(iwave == 1)THEN
392C add afterflow velocity, no afterflow velocity for plane wave
393 IF(ipres == 1) THEN
394 DO i=1,nel
395 ratio = dsc/dist(i)
396 pmax = rflow(6)*ratio**apmax
397 theta = rflow(7)*ratio**atheta
398 vi(i) = vi(i) + (pmax-pp(i))*cosg(i)*theta/(rho*dist(i))
399 ENDDO
400 ELSEIF(ipres == 2) THEN
401 DO i=1,nel
402 pti(i) = pti(i) + pp(i)*dt1
403 vi(i) = vi(i) + pti(i)*cosg(i)/(rho*dist(i))
404 ENDDO
405 ENDIF
406 ENDIF
407 ENDIF
408C-----------------------------------------------
409C 2. Pression Onde Reflechie
410C-----------------------------------------------
411 IF(freesurf == 2) THEN
412 DO i=1,nel
413 off(i) = zero
414 tta = tt-ta(i+nel)
415 IF(tta > zero) off(i) = one
416 ENDDO
417 IF(ipres == 1) THEN
418 DO i=1,nel
419 j = i+nel
420 tta = tt-ta(j)
421 ratio = dsc/dist(j)
422 pmax = rflow(6)*ratio**apmax
423 theta = rflow(7)*ratio**atheta
424 pr(i) = -pmax*exp(-tta/theta)*off(i)
425 dprdt(i) = -pr(i)*cosg(j)/theta
426 ENDDO
427 ELSEIF(ipres == 2) THEN
428 DO i=1,nel
429 j = i+nel
430 tta = tt-ta(j)
431 ratio = dsc/dist(j)
432 IF(ifpres > 0) THEN
433 ismooth = npc(2*nfunct+ifpres+1)
434!! PR(I) = -SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
435 IF (ismooth == 0) THEN
436 pr(i) = -sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
437 ELSE IF (ismooth > 0) THEN
438 pr(i) = -sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
439 ELSE
440 ismooth = -ismooth ! function id
441 CALL python_call_funct1d(python, ismooth,tta,pr(i))
442 pr(i) = -pr(i)*sfpres*ratio*off(i)
443 CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
444 ENDIF ! IF (ISMOOTH == 0)
445 dprdt(i) = dppdt*cosg(j)*off(i)
446 ELSE
447 pr(i) = -sfpres*ratio*off(i)
448 dprdt(i) = zero
449 ENDIF
450 ENDDO
451 ENDIF
452C-----------------------------------------------
453C Vitesse onde reflechie
454C-----------------------------------------------
455 DO i=1,nel
456 vr(i)=pr(i)*cosg(i+nel)/rhoc
457 ENDDO
458 ELSE ! pas de surface libre
459 DO i=1,nel
460 pr(i) = zero
461 vr(i) = zero
462 dprdt(i) = zero
463 ENDDO
464 ENDIF
465C-----------------------------------------------
466C 3. Scattered pressure Schema Euler ordre 1
467C-----------------------------------------------
468 IF(kform == 1) THEN
469C DAA-1
470 IF(nel < nelmax) THEN
471 DO i=1,nel
472 dpmdt(i) = rhoc*accf(i)
473 vect(i) = rhoc*areaf(i)*(pm(i)-rhoc*(vi(i)+vr(i)))
474 ENDDO
475 vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
476 DO i=1,nel
477 dpmdt(i) = dpmdt(i) - vect(i)
478 ENDDO
479 IF(integr == 1) THEN
480 DO i=1,nel
481 pm(i) = pm(i) + dpmdt(i)*dt1
482 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
483 ENDDO
484
485 ELSEIF(integr == 2) THEN
486C Prediction
487 DO i=1,nel
488 ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
489 ENDDO
490C Correction
491 DO i=1,nel
492 dpmdt(i) = rhoc*accf(i)
493 vect(i) = rhoc*areaf(i)*ps(i)
494 ENDDO
495 vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
496 DO i=1,nel
497 dpmdt(i) = dpmdt(i) - vect(i)
498 ENDDO
499 DO i=1,nel
500 pm(i) = pm(i) + dpmdt(i)*dt1
501 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
502 ENDDO
503 ENDIF
504C
505 ELSE ! resolution systeme MFLE*X=Y
506C Initialisation de la process grid
507 CALL sl_init(ictxt, nprow, npcol)
508 CALL blacs_gridinfo(ictxt, nprow, npcol, myrow, mycol)
509C MXLLDM = NUMROC(NEL, NBLOC, MYROW, 1, NPROW)
510 mxlldm = nel_loc
511 CALL descinit(desca, nel, nel, nbloc, nbloc, 0, 0, ictxt, mxlldm, info)
512 CALL descinit(descb, nel, 1, nbloc, nbloc, 0, 0, ictxt, mxlldm, info)
513
514 IF(dt1 == zero) THEN
515C Compute local sub-matrix
516 k=0
517 DO i=1,nel
518 IF(ibufelr(i) == 0) cycle
519 k=k+1
520 l=0
521 DO j=1,nel
522 IF(ibufelc(j) == 0) cycle
523 l=l+1
524 mfle_l(k,l) = mfle(i,j)
525 ENDDO
526 ENDDO
527C L U decomposition
528 CALL pdgetrf(nel, nel, mfle_l, 1, 1, desca, ipiv_l, info)
529 ENDIF
530
531 DO i=1,nel
532 vect(i) = two*ssp*(rhoc*(vi(i)+vr(i))-pm(i))
533 ENDDO
534 CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
535
536 vect(1:nel)=zero
537 IF(ibufelc(1) > 0) THEN
538 k=0
539 DO i=1,nel
540 IF(ibufelr(i) == 0) cycle
541 k=k+1
542 vect(k) = vect1(i)
543 ENDDO
544 ENDIF
545
546C Resolution
547 CALL pdgetrs('N', nel, 1, mfle_l, 1, 1, desca, ipiv_l, vect, 1, 1, descb, info)
548
549C Echange
550 ALLOCATE(sbuf(nel), rbuf(nel))
551 sbuf(1:nel)=zero
552 rbuf(1:nel)=zero
553 IF(ibufelc(1) > 0) THEN
554 k=0
555 DO i=1,nel
556 IF(ibufelr(i) == 0) cycle
557 k=k+1
558 sbuf(i) = vect(k)
559 ENDDO
560 ENDIF
561 CALL spmd_fl_sum(sbuf, nel, rbuf)
562 DO i=1,nel
563 vect(i)=rbuf(i)
564 ENDDO
565
566 vect1(1:nel) = matmul(cbem(1:nel,1:nel), vect(1:nel))
567 DO i=1,nel
568 IF(pp(i) /= zero) THEN
569 dpmdt(i) = rhoc*accf(i) + vect(i)
570 ELSE
571 dpmdt(i)=zero
572 ENDIF
573 ENDDO
574
575 IF(integr == 1) THEN
576 DEALLOCATE(sbuf, rbuf)
577 DO i=1,nel
578 pm(i) = pm(i) + dpmdt(i)*dt1
579 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
580 ENDDO
581
582 ELSEIF(integr == 2) THEN
583C Prediction
584 DO i=1,nel
585 ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
586 ENDDO
587C Correction
588 DO i=1,nel
589 vect(i) = -two*ssp*ps(i)
590 ENDDO
591 CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
592 vect(1:nel)=zero
593 IF(ibufelc(1) > 0) THEN
594 k=0
595 DO i=1,nel
596 IF(ibufelr(i) == 0) cycle
597 k=k+1
598 vect(k) = vect1(i)
599 ENDDO
600 ENDIF
601
602C Resolution
603 CALL pdgetrs('N', nel, 1, mfle_l, 1, 1, desca, ipiv_l, vect, 1, 1, descb, info)
604C Echange
605 sbuf(1:nel)=zero
606 rbuf(1:nel)=zero
607 IF(ibufelc(1) > 0) THEN
608 k=0
609 DO i=1,nel
610 IF(ibufelr(i) == 0) cycle
611 k=k+1
612 sbuf(i) = vect(k)
613 ENDDO
614 ENDIF
615 CALL spmd_fl_sum(sbuf, nel, rbuf)
616 DO i=1,nel
617 vect(i)=rbuf(i)
618 ENDDO
619 DEALLOCATE(sbuf, rbuf)
620
621 vect1(1:nel) = matmul(cbem(1:nel,1:nel), vect(1:nel))
622 DO i=1,nel
623 dpmdt(i) = rhoc*accf(i) + vect1(i)
624 ENDDO
625 DO i=1,nel
626 pm(i) = pm(i) + dpmdt(i)*dt1
627 ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
628 ENDDO
629 ENDIF
630
631 CALL blacs_gridexit(ictxt)
632 ENDIF
633C
634 ELSEIF(kform == 2) THEN
635C High Frequency Approximation
636 DO i=1,nel
637 dpmdt(i) = rhoc*accf(i)
638 ENDDO
639 DO i=1,nel
640 pm(i) = pm(i) + dpmdt(i)*dt1
641 ps(i) = pm(i) - rhoc*vi(i)
642 ENDDO
643C
644 ELSEIF(kform == 3) THEN
645C Virtual Mass Approximation
646 DO i=1,nel
647 ps(i) = zero
648 DO j=1,nel
649 ps(i) = ps(i) + mfle(i,j)*(accf(j)-(dpidt(j)+dprdt(j))/rhoc)
650 ENDDO
651 ps(i) = ps(i)/areaf(i)
652 ENDDO
653 ENDIF
654C-----------------------------------------------
655C 4. gravity
656C-----------------------------------------------
657 gravity = zero
658 IF(grav_id > 0) THEN
659 fcy = agrv(1,grav_id)
660 fcx = agrv(2,grav_id)
661 ifunc = igrv(3,grav_id)
662 isens = 0
663 DO k=1,nsensor
664 IF(igrv(6,grav_id) == sensor_tab(k)%SENS_ID) isens=k ! do it in starter !!!
665 ENDDO
666 IF(isens==0)THEN
667 ts = tt
668 ELSE
669 ts = tt-sensor_tab(isens)%TSTART
670 ENDIF
671 ismooth = 0
672 IF (ifunc > 0) THEN
673 ismooth = npc(2*nfunct+ifunc+1)
674!! GRAVITY = FCY*FINTER(IFUNC,TS*FCX,NPC,TF,DYDX)
675 IF (ismooth == 0) THEN
676 gravity = fcy*finter(ifunc,ts*fcx,npc,tf,dydx)
677 ELSE IF(ismooth > 0) THEN
678 gravity = fcy*finter_smooth(ifunc,ts*fcx,npc,tf,dydx)
679 ELSE
680 ismooth = -ismooth ! function id
681 CALL python_call_funct1d(python, ismooth,tta,gravity)
682 gravity = gravity*fcy
683 ENDIF ! IF (ISMOOTH == 0)
684 ELSE
685 gravity = fcy
686 ENDIF
687 ENDIF
688C-----------------------------------------------
689C 5. compute total pressure forces
690C-----------------------------------------------
691 DO i=1,nel
692 ptot = pp(i)+ps(i)+pr(i)+rho*hh(i)*abs(gravity)
693 fel(i) = areaf(i)*max(ptot,pmin)
694 ENDDO
695C
696 IF(jform == 1) THEN
697 nel_l=0
698 DO i=1,nel
699 i1 = elem(1,i)
700 i2 = elem(2,i)
701 i3 = elem(3,i)
702 n1 = ibuf(i1)
703 n2 = ibuf(i2)
704 n3 = ibuf(i3)
705 IF(n1/=0.OR.n2/=0.OR.n3/=0) THEN
706 nel_l = nel_l+1
707 nebuf(nel_l) = i
708 ENDIF
709 ENDDO
710 DO k=1,nel_l
711 i = nebuf(k)
712 i1 = elem(1,i)
713 i2 = elem(2,i)
714 i3 = elem(3,i)
715 n1 = ibuf(i1)
716 n2 = ibuf(i2)
717 n3 = ibuf(i3)
718 nrx= normal(1,i)
719 nry= normal(2,i)
720 nrz= normal(3,i)
721 ff = fac1*third*fel(i)
722 IF(n1/=0) THEN
723 a(1,n1) = a(1,n1) + ff*nrx
724 a(2,n1) = a(2,n1) + ff*nry
725 a(3,n1) = a(3,n1) + ff*nrz
726 wfext = wfext+ff*dt1*(nrx*v(1,n1)+nry*v(2,n1)+nrz*v(3,n1))/cnp(i1)
727 ENDIF
728 IF(n2/=0) THEN
729 a(1,n2) = a(1,n2) + ff*nrx
730 a(2,n2) = a(2,n2) + ff*nry
731 a(3,n2) = a(3,n2) + ff*nrz
732 wfext = wfext+ff*dt1*(nrx*v(1,n2)+nry*v(2,n2)+nrz*v(3,n2))/cnp(i2)
733 ENDIF
734 IF(n3/=0) THEN
735 a(1,n3) = a(1,n3) + ff*nrx
736 a(2,n3) = a(2,n3) + ff*nry
737 a(3,n3) = a(3,n3) + ff*nrz
738 wfext = wfext+ff*dt1*(nrx*v(1,n3)+nry*v(2,n3)+nrz*v(3,n3))/cnp(i3)
739 ENDIF
740 ENDDO
741 ELSEIF(jform == 2) THEN
742 nel_l=0
743 DO i=1,nel
744 i1 = elem(1,i)
745 i2 = elem(2,i)
746 i3 = elem(3,i)
747 i4 = elem(4,i)
748 n1 = ibuf(i1)
749 n2 = ibuf(i2)
750 n3 = ibuf(i3)
751 n4 = ibuf(i4)
752 IF(n1/=0.OR.n2/=0.OR.n3/=0.OR.n4/=0) THEN
753 nel_l = nel_l+1
754 nebuf(nel_l) = i
755 ENDIF
756 ENDDO
757 wf(1)=fac1*fourth
758 wf(2)=fac1*third
759 DO k=1,nel_l
760 i = nebuf(k)
761 i1 = elem(1,i)
762 i2 = elem(2,i)
763 i3 = elem(3,i)
764 i4 = elem(4,i)
765 i5 = elem(5,i)
766 n1 = ibuf(i1)
767 n2 = ibuf(i2)
768 n3 = ibuf(i3)
769 n4 = ibuf(i4)
770 nrx= normal(1,i)
771 nry= normal(2,i)
772 nrz= normal(3,i)
773 ff = wf(i5)*fel(i)
774 IF(n1/=0) THEN
775 a(1,n1) = a(1,n1) + ff*nrx
776 a(2,n1) = a(2,n1) + ff*nry
777 a(3,n1) = a(3,n1) + ff*nrz
778 wfext = wfext+ff*dt1*(nrx*v(1,n1)+nry*v(2,n1)+nrz*v(3,n1))/cnp(i1)
779 ENDIF
780 IF(n2/=0) THEN
781 a(1,n2) = a(1,n2) + ff*nrx
782 a(2,n2) = a(2,n2) + ff*nry
783 a(3,n2) = a(3,n2) + ff*nrz
784 wfext = wfext+ff*dt1*(nrx*v(1,n2)+nry*v(2,n2)+nrz*v(3,n2))/cnp(i2)
785 ENDIF
786 IF(n3/=0) THEN
787 a(1,n3) = a(1,n3) + ff*nrx
788 a(2,n3) = a(2,n3) + ff*nry
789 a(3,n3) = a(3,n3) + ff*nrz
790 wfext = wfext+ff*dt1*(nrx*v(1,n3)+nry*v(2,n3)+nrz*v(3,n3))/cnp(i3)
791 ENDIF
792 IF(n4/=0.AND.i5==1) THEN
793 a(1,n4) = a(1,n4) + ff*nrx
794 a(2,n4) = a(2,n4) + ff*nry
795 a(3,n4) = a(3,n4) + ff*nrz
796 wfext = wfext+ff*dt1*(nrx*v(1,n4)+nry*v(2,n4)+nrz*v(3,n4))/cnp(i4)
797 ENDIF
798 ENDDO
799 ENDIF
800C-----------------------------------------------------------
801C 5. Gauge
802C-----------------------------------------------------------
803 IF(jform == 2) THEN
804 DO i=1,nbgauge
805 gauge(30,i)=zero
806 IF(lgauge(1,i) == 0 .AND. lgauge(3,i) < 0 ) THEN
807 i1=shell_ga(i)
808 IF(i1 > 0) THEN
809 ptot=pp(i1)+ps(i1)+pr(i1)+rho*hh(i1)*abs(gravity)
810 gauge(30,i)=max(ptot,pmin)
811 ENDIF
812 ENDIF
813 ENDDO
814 ENDIF
815
816#endif
817 RETURN
sl_init
subroutine sl_init(ictxt, nprow, npcol)
Definition SL_init.f:2
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
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
pdgetrs
subroutine pdgetrs(trans, n, nrhs, a, ia, ja, desca, ipiv, b, ib, jb, descb, info)
Definition mpi.f:1171
pdgetrf
subroutine pdgetrf(m, n, a, ia, ja, desca, ipiv, info)
Definition mpi.f:914
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition mpi.f:777
blacs_gridexit
subroutine blacs_gridexit(cntxt)
Definition mpi.f:762
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition mpi.f:786
spmd_fl_sum
subroutine spmd_fl_sum(lsum, len, lsumt)
Definition spmd_fl_sum.F:37