OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sponfv.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!|| sponfv ../engine/source/elements/sph/sponfv.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| spmd_sphgeta ../engine/source/mpi/elements/spmd_sph.F
29!|| spmd_sphgetd ../engine/source/mpi/elements/spmd_sph.F
30!|| weight0 ../engine/source/elements/sph/weight.F
31!|| weight1 ../engine/source/elements/sph/weight.F
32!||--- uses -----------------------------------------------------
33!|| sph_work_mod ../common_source/modules/mat_elem/sph_work.F90
34!|| sphbox ../engine/share/modules/sphbox.F
35!||====================================================================
36 SUBROUTINE sponfv(X ,V ,A ,D ,MS ,
37 2 SPBUF ,ITAB ,KXSP ,IXSP ,NOD2SP ,
38 3 NPC ,PLD ,ISPHIO ,VSPHIO ,IPART ,
39 4 IPARTSP ,WASPACT ,WA ,VNORMAL ,SPH_WORK,
40 5 WFEXT )
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE sphbox
45 USE sph_work_mod
46C-----------------------------------------------
47C I m p l i c i t T y p e s
48C-----------------------------------------------
49#include "implicit_f.inc"
50#include "comlock.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "com01_c.inc"
55#include "com06_c.inc"
56#include "com08_c.inc"
57#include "sphcom.inc"
58#include "scr17_c.inc"
59C-----------------------------------------------
60C D u m m y A r g u m e n t s
61C-----------------------------------------------
62 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),NPC(*),
63 . ISPHIO(NISPHIO,*),IPART(LIPART1,*),IPARTSP(*),
64 . WASPACT(*)
65 my_real X(3,*) ,V(3,*) ,A(3,*) ,D(3,*) ,MS(*) ,SPBUF(NSPBUF,*) ,PLD(*) ,VSPHIO(*), WA(*), VNORMAL(3,*)
66 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
67C-----------------------------------------------
68C L o c a l V a r i a b l e s
69C-----------------------------------------------
70 INTEGER I, ITYPE,
71 . II,IPT,JJ,NPF,IFVITS,
72 . NS,N,INOD,IACTIVE,
73 . IPPV,J,M,JNOD,IMPOSE,JMPOSE,
74 . nvois,ij,np,k,jmpose2,nn
75 my_real
76 . wfextt,t05,
77 . pentv,vx,vy,vz,vn,vt,ux,uy,uz,un1,nx,ny,nz,
78 . ps,
79 . xi,yi,zi,xj,yj,zj,dmin,dd,dtinv,
80 . di,rhoi,dj,rhoj,dij,
81 . vxi,vyi,vzi,vxj,vyj,vzj,
82 . vj,vjx,vjy,vjz,
83 . wght,wgrad(3),wgrdx,wgrdy,wgrdz,
84 . dxx,dxy,dxz,dyx,dyy,dyz,dzx,dzy,dzz,dt1d2,
85 . exx,exy,exz,eyx,eyy,eyz,ezx,ezy,ezz,
86 . alphai,alphaxi,alphayi,alphazi,alphai2,xp,yp,zp
87
88 LOGICAL lBOOL
89
90 my_real ,DIMENSION(:,:), ALLOCATABLE :: DSPHR
91 TYPE(SPH_WORK_) :: SPH_WORK
92C-----------------------------------------------
93C inlets.
94C-----------------------------------------------
95 t05=tt + half*dt2
96 xp = 0
97 yp = 0
98 zp = 0
99
100 DO i=1,nsphio
101 itype=isphio(1,i)
102 IF(itype==1)THEN
103C------
104 ifvits=isphio(8,i)
105 vn =zero
106 IF(ifvits/=0)THEN
107 npf = (npc(ifvits+1)-npc(ifvits))/2
108 ii = npc(ifvits)
109 IF (t05<=pld(ii)) THEN
110 pentv=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
111 vn =pld(ii+1)+pentv*(t05-pld(ii))
112 ELSEIF (t05>=pld(ii+2*(npf-1))) THEN
113 jj=ii+2*(npf-1)
114 pentv=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
115 vn =pld(jj+1)+max(-pld(jj+1),pentv*(t05-pld(jj)))
116 ELSE
117 DO ipt=1,npf-1
118 IF (pld(ii)<=t05.AND.t05<=pld(ii+2)) THEN
119 pentv=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
120 vn =pld(ii+1)+pentv*(t05-pld(ii))
121 GOTO 10
122 ENDIF
123 ii=ii+2
124 ENDDO
125 10 CONTINUE
126 ENDIF
127 ENDIF
128 wa(i)=vn
129 ENDIF
130 ENDDO
131
132C-----
133 wfextt=zero
134 dtinv =one/dt12
135C
136 DO ns=1,nsphact
137 n=waspact(ns)
138 impose=kxsp(2,n)/(ngroup+1)
139 IF(impose/=0)THEN
140 itype=isphio(1,impose)
141 IF(itype==1)THEN
142C------
143 vn=wa(impose)
144C------
145 inod=kxsp(3,n)
146 ux=v(1,inod)+a(1,inod)*dt12
147 uy=v(2,inod)+a(2,inod)*dt12
148 uz=v(3,inod)+a(3,inod)*dt12
149 nx=vnormal(1,n)
150 ny=vnormal(2,n)
151 nz=vnormal(3,n)
152 un1=ux*nx+uy*ny+uz*nz
153 vx=ux+(vn-un1)*nx
154 vy=uy+(vn-un1)*ny
155 vz=uz+(vn-un1)*nz
156 wfextt=wfextt+half*ms(inod)*
157 . ((vx*vx+vy*vy+vz*vz)-(ux*ux+uy*uy+uz*uz))
158 a(1,inod)=(vx-v(1,inod))*dtinv
159 a(2,inod)=(vy-v(2,inod))*dtinv
160 a(3,inod)=(vz-v(3,inod))*dtinv
161 ENDIF
162 ENDIF
163 ENDDO
164
165C-----------------------------------------------
166C Comm D et A sur cellules remotes
167C-----------------------------------------------
168 IF(nspmd>1)THEN
169 ALLOCATE(sph_work%ASPHR(3,nsphr))
170 ALLOCATE(dsphr(12,nsphr))
171 CALL spmd_sphgeta(kxsp,spbuf,a,sph_work%ASPHR)
172 CALL spmd_sphgetd(kxsp,ixsp,isphio,x,waspact,nod2sp,
173 . spbuf,v,a,sph_work%ASPHR,dsphr)
174 END IF
175C-------------------------------------------
176C general outlet & silent boundary.
177C-------------------------------------------
178
179 DO ns=1,nsphact
180 n=waspact(ns)
181 impose=kxsp(2,n)/(ngroup+1)
182 lbool=.false.
183 IF(impose /= 0)THEN
184 IF(isphio(1,impose)==2.OR.isphio(1,impose)==3)lbool=.true.
185 ENDIF
186 IF(lbool)THEN
187 inod=kxsp(3,n)
188 xi=x(1,inod)
189 yi=x(2,inod)
190 zi=x(3,inod)
191C-------
192C plus proche voisin en amont de l'outlet => IPPV.
193 ippv=0
194 dmin=1.e+20
195 DO k=1,kxsp(4,n)
196 jnod=ixsp(k,n)
197
198 IF(jnod>0)THEN
199 m =nod2sp(jnod)
200 jmpose=kxsp(2,m)/(ngroup+1)
201 lbool=.false.
202 IF(jmpose == 0)THEN
203 lbool=.true.
204 ELSE
205 IF(isphio(1,jmpose) == 1)lbool=.true.
206 ENDIF
207 IF(lbool)THEN
208 xj =x(1,jnod)
209 yj =x(2,jnod)
210 zj =x(3,jnod)
211 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
212 IF(dd<dmin)THEN
213 ippv=jnod
214 dmin=dd
215 ENDIF
216 ENDIF
217 ELSE
218 nn = -jnod
219 jmpose = nint(xsphr(12,nn))
220 IF(jmpose>0)THEN
221 jmpose2=isphio(1,jmpose)
222 ELSE
223 jmpose2=0
224 ENDIF
225 IF(jmpose2==0.OR.jmpose2==1)THEN
226 xj =xsphr(3,nn)
227 yj =xsphr(4,nn)
228 zj =xsphr(5,nn)
229 dd =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
230 IF(dd<dmin)THEN
231 ippv=jnod
232 dmin=dd
233 ENDIF
234 ENDIF
235 ENDIF
236 ENDDO
237C-------
238C calcule grad.v en IPPV (approximation sph corrigee en amont de IPPV).
239 IF(ippv>0)THEN
240 np=nod2sp(ippv)
241 xp=x(1,ippv)
242 yp=x(2,ippv)
243 zp=x(3,ippv)
244 di =spbuf(1,np)
245 rhoi=spbuf(2,np)
246
247 CALL weight0(xp,yp,zp,xp,yp,zp,di,wght)
248 vj=spbuf(12,np)/max(em20,rhoi)
249 alphai=vj*wght
250 alphaxi=zero
251 alphayi=zero
252 alphazi=zero
253
254 DO j=1,kxsp(4,np)
255 jnod=ixsp(j,np)
256
257c partie JNOD > 0
258 IF(jnod>0)THEN ! particule locale
259 m=nod2sp(jnod)
260 jmpose=kxsp(2,m)/(ngroup+1)
261 lbool=.false.
262 IF(jmpose == 0)THEN
263 lbool=.true.
264 ELSE
265 IF(isphio(1,jmpose) == 1)lbool=.true.
266 ENDIF
267 IF(lbool)THEN
268 dj =spbuf(1,m)
269 xj =x(1,jnod)
270 yj =x(2,jnod)
271 zj =x(3,jnod)
272 dij =(dj+di)*half
273 rhoj=spbuf(2,m)
274 CALL weight1(xp,yp,zp,xj,yj,zj,dij,wght,wgrad)
275 vj=spbuf(12,m)/max(em20,rhoj)
276 alphai =alphai +vj*wght
277 alphaxi=alphaxi+vj*wgrad(1)
278 alphayi=alphayi+vj*wgrad(2)
279 alphazi=alphazi+vj*wgrad(3)
280 ENDIF
281 ELSE ! particule remote
282 nn = -jnod
283 jmpose = nint(xsphr(12,nn))
284
285 IF(jmpose>0)THEN
286 jmpose2=isphio(1,jmpose)
287 ELSE
288 jmpose2=0
289 ENDIF
290 IF(jmpose2==0.OR.jmpose2==1)THEN
291 dj =xsphr(2,nn)
292 xj =xsphr(3,nn)
293 yj =xsphr(4,nn)
294 zj =xsphr(5,nn)
295 dij =(dj+di)*half
296 rhoj=xsphr(7,nn)
297 CALL weight1(xp,yp,zp,xj,yj,zj,dij,wght,wgrad)
298 vj=xsphr(8,nn)/max(em20,rhoj)
299 alphai =alphai +vj*wght
300 alphaxi=alphaxi+vj*wgrad(1)
301 alphayi=alphayi+vj*wgrad(2)
302 alphazi=alphazi+vj*wgrad(3)
303 ENDIF
304 ENDIF
305 ENDDO
306C------
307 alphai =one/max(em20,alphai)
308 alphai2=alphai*alphai
309 alphaxi=-alphaxi*alphai2
310 alphayi=-alphayi*alphai2
311 alphazi=-alphazi*alphai2
312C------
313 vx =v(1,ippv)+dt12*a(1,ippv)
314 vy =v(2,ippv)+dt12*a(2,ippv)
315 vz =v(3,ippv)+dt12*a(3,ippv)
316 dxx=zero
317 dxy=zero
318 dxz=zero
319 dyx=zero
320 dyy=zero
321 dyz=zero
322 dzx=zero
323 dzy=zero
324 dzz=zero
325
326 DO j=1,kxsp(4,np)
327 jnod=ixsp(j,np)
328 IF(jnod>0)THEN
329 m=nod2sp(jnod)
330 jmpose=kxsp(2,m)/(ngroup+1)
331 lbool=.false.
332 IF(jmpose == 0)THEN
333 lbool=.true.
334 ELSE
335 IF(isphio(1,jmpose) == 1)lbool=.true.
336 ENDIF
337 IF(lbool)THEN
338 dj =spbuf(1,m)
339 xj =x(1,jnod)
340 yj =x(2,jnod)
341 zj =x(3,jnod)
342 dij =(dj+di)*half
343 rhoj=spbuf(2,m)
344 CALL weight1(xp,yp,zp,xj,yj,zj,dij,wght,wgrad)
345 wgrdx=wgrad(1)*alphai+wght*alphaxi
346 wgrdy=wgrad(2)*alphai+wght*alphayi
347 wgrdz=wgrad(3)*alphai+wght*alphazi
348 vj=spbuf(12,m)/max(em20,rhoj)
349 vxj =v(1,jnod)+dt12*a(1,jnod)
350 vyj =v(2,jnod)+dt12*a(2,jnod)
351 vzj =v(3,jnod)+dt12*a(3,jnod)
352 vjx=vj*(vxj-vx)
353 vjy=vj*(vyj-vy)
354 vjz=vj*(vzj-vz)
355 dxx=dxx+vjx*wgrdx
356 dxy=dxy+vjx*wgrdy
357 dxz=dxz+vjx*wgrdz
358 dyx=dyx+vjy*wgrdx
359 dyy=dyy+vjy*wgrdy
360 dyz=dyz+vjy*wgrdz
361 dzx=dzx+vjz*wgrdx
362 dzy=dzy+vjz*wgrdy
363 dzz=dzz+vjz*wgrdz
364 ENDIF
365 ELSE
366 nn=-jnod
367 jmpose = nint(xsphr(12,nn))
368 IF(jmpose>0)THEN
369 jmpose2=isphio(1,jmpose)
370 ELSE
371 jmpose2=0
372 ENDIF
373 IF(jmpose2==0.OR.jmpose2==1)THEN
374 dj =xsphr(2,nn)
375 xj =xsphr(3,nn)
376 yj =xsphr(4,nn)
377 zj =xsphr(5,nn)
378 dij =(dj+di)*half
379 rhoj=xsphr(7,nn)
380 CALL weight1(xp,yp,zp,xj,yj,zj,dij,wght,wgrad)
381 wgrdx=wgrad(1)*alphai+wght*alphaxi
382 wgrdy=wgrad(2)*alphai+wght*alphayi
383 wgrdz=wgrad(3)*alphai+wght*alphazi
384 vj=xsphr(8,nn)/max(em20,rhoj)
385 vxj =xsphr(9,nn)+dt12*sph_work%ASPHR(1,nn)
386 vyj =xsphr(10,nn)+dt12*sph_work%ASPHR(2,nn)
387 vzj =xsphr(11,nn)+dt12*sph_work%ASPHR(3,nn)
388 vjx=vj*(vxj-vx)
389 vjy=vj*(vyj-vy)
390 vjz=vj*(vzj-vz)
391 dxx=dxx+vjx*wgrdx
392 dxy=dxy+vjx*wgrdy
393 dxz=dxz+vjx*wgrdz
394 dyx=dyx+vjy*wgrdx
395 dyy=dyy+vjy*wgrdy
396 dyz=dyz+vjy*wgrdz
397 dzx=dzx+vjz*wgrdx
398 dzy=dzy+vjz*wgrdy
399 dzz=dzz+vjz*wgrdz
400 ENDIF
401 ENDIF
402 ENDDO
403C------
404 ELSEIF(ippv<0)then! fin traitement IPPV > 0
405c traitement IPPV negatif on utilise infos recuperres ds routine comm
406 dxx = dsphr(1,-ippv)
407 dxy = dsphr(2,-ippv)
408 dxz = dsphr(3,-ippv)
409 dyx = dsphr(4,-ippv)
410 dyy = dsphr(5,-ippv)
411 dyz = dsphr(6,-ippv)
412 dzx = dsphr(7,-ippv)
413 dzy = dsphr(8,-ippv)
414 dzz = dsphr(9,-ippv)
415 vx = dsphr(10,-ippv)
416 vy = dsphr(11,-ippv)
417 vz = dsphr(12,-ippv)
418 xp=xsphr(3,-ippv)
419 yp=xsphr(4,-ippv)
420 zp=xsphr(5,-ippv)
421c ELSE ! cas IPPV=0 error
422c print*,'IPPV nul ERROR'
423 endif! fin traitement IPPV < 0
424
425 vx=vx+(dxx*(xi-xp)+dxy*(yi-yp)+dxz*(zi-zp))
426 vy=vy+(dyx*(xi-xp)+dyy*(yi-yp)+dyz*(zi-zp))
427 vz=vz+(dzx*(xi-xp)+dzy*(yi-yp)+dzz*(zi-zp))
428
429 ps=vx*vnormal(1,n)+vy*vnormal(2,n)+vz*vnormal(3,n)
430
431 IF(ps<0.)THEN
432C impose une vitesse sortante uniquement.
433 vx=vx-ps*vnormal(1,n)
434 vy=vy-ps*vnormal(2,n)
435 vz=vz-ps*vnormal(3,n)
436 ENDIF
437 ux=v(1,inod)+a(1,inod)*dt12
438 uy=v(2,inod)+a(2,inod)*dt12
439 uz=v(3,inod)+a(3,inod)*dt12
440
441 vt=vx*vx+vy*vy+vz*vz
442 wfextt=wfextt+half*ms(inod)*(vt-(ux*ux+uy*uy+uz*uz))
443 a(1,inod)=(vx-v(1,inod))*dtinv
444 a(2,inod)=(vy-v(2,inod))*dtinv
445 a(3,inod)=(vz-v(3,inod))*dtinv
446 ENDIF
447 ENDDO
448C-------------------------------------------
449!$OMP ATOMIC
450 wfext=wfext+wfextt
451C-------------------------------------------
452 IF(nspmd>1)THEN
453 DEALLOCATE(sph_work%ASPHR,dsphr)
454 END IF
455C-------------------------------------------
456 RETURN
457 END
max
#define max(a, b)
Definition macros.h:21
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
spmd_sphgetd
subroutine spmd_sphgetd(kxsp, ixsp, isphio, x, waspact, nod2sp, spbuf, v, a, asphr, dsphr)
Definition spmd_sph.F:1929
spmd_sphgeta
subroutine spmd_sphgeta(kxsp, spbuf, a, asphr)
Definition spmd_sph.F:822
sponfv
subroutine sponfv(x, v, a, d, ms, spbuf, itab, kxsp, ixsp, nod2sp, npc, pld, isphio, vsphio, ipart, ipartsp, waspact, wa, vnormal, sph_work, wfext)
Definition sponfv.F:41
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79
weight0
subroutine weight0(xi, yi, zi, xj, yj, zj, h, w)
Definition weight.F:34