OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spforcp.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!|| spforcp ../engine/source/elements/sph/spforcp.F
25!||--- called by ------------------------------------------------------
26!|| forintp ../engine/source/elements/forintp.F
27!||--- calls -----------------------------------------------------
28!|| weight1 ../engine/source/elements/sph/weight.F
29!||--- uses -----------------------------------------------------
30!|| sphbox ../engine/share/modules/sphbox.F
31!||====================================================================
32 SUBROUTINE spforcp(
33 1 PM ,GEO ,X ,V ,MS ,
34 2 SPBUF ,ITAB ,PLD ,BUFMAT ,BUFGEO ,
35 3 PARTSAV ,FSAV ,DT2T ,IPARG ,NPC ,
36 4 KXSP ,IXSP ,NOD2SP ,NELTST ,ITYPTST ,
37 5 IPART ,IPARTSP ,ISPCOND ,XFRAME ,ISPSYM ,
38 6 XSPSYM ,VSPSYM ,WA ,WASIGSM ,WACOMP ,
39 7 WSMCOMP ,WASPACT ,WAR ,STAB ,WFEXT)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE sphbox
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "vect01_c.inc"
53#include "com06_c.inc"
54#include "com08_c.inc"
55#include "sphcom.inc"
56#include "param_c.inc"
57#include "scr02_c.inc"
58#include "scr17_c.inc"
59#include "task_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
64 . IPART(LIPART1,*) ,IPARTSP(*), NPC(*), IPARG(NPARG,*),
65 . NELTST,ITYPTST,ITAB(*) ,ISPCOND(NISPCOND,*),
66 . ISPSYM(NSPCOND,*), WASPACT(*)
67 my_real
68 . X(3,*) ,V(3,*) ,MS(*) ,
69 . PM(NPROPM,*), GEO(NPROPG,*),
70 . bufmat(*) ,bufgeo(*) ,pld(*) ,fsav(nthvki,*) ,
71 . spbuf(nspbuf,*) ,partsav(npsav,*) ,dt2t ,
72 . xframe(nxframe,*) ,xspsym(3,*) ,vspsym(3,*), wa(kwasph,*),
73 . wasigsm(6,*), wacomp(16,*), wsmcomp(6,*), war(10,*), stab(7,*)
74 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER N, INOD, JNOD, J, NVOIS, M, IPRT, IPROP, I,
79 . NVOISS,NC,SM,JS,NN,KS,NR
80 my_real
81 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,
82 . sxx,sxy,sxz,syy,syz,szz,
83 . txx,txy,txz,tyy,tyz,tzz,
84 . ax,ay,az,bx,by,bz,fx,fy,fz,
85 . wght,wgrad(3),
86 . dij,mm,
87 . vxi,vyi,vzi,vxj,vyj,vzj,muij,muij2,pij,ssp,
88 . fact,qa,qb,
89 . divvi,divvj,rotvi,rotvj,fi,fj,
90 . fvx,fvy,fvz,wvis,
91 . stij,sfac,dldt,l,cij,dzeta,wnorm,
92 . vi,vj,vij,
93 . ssp2i,ssp2j,stii,
94 . alphai,alphaxi,alphayi,alphazi,
95 . betaxxi,betayxi,betazxi,
96 . betaxyi,betayyi,betazyi,
97 . betaxzi,betayzi,betazzi,
98 . alphaj,alphaxj,alphayj,alphazj,
99 . betaxj,betayj,betazj,
100 . betaxxj,betayxj,betazxj,
101 . betaxyj,betayyj,betazyj,
102 . betaxzj,betayzj,betazzj,
103 . wgrdx,wgrdy,wgrdz,wgrd(3),drhojdr,drhoidr,
104 . voli,volj,dxij,nxij,
105 . cx,cy,cz,dx,dy,dz,tx,ty,tz,wfextt, ww, wi, wr
106C-----------------------------------------------
107 wfextt=zero
108 DO 10 i=lft,llt
109 n =nft+i
110 IF(kxsp(2,n)<=0.AND.isph2sol==0)GOTO 10
111 IF(kxsp(2,n)==0.AND.isph2sol/=0)GOTO 10
112 inod =kxsp(3,n)
113 xi=x(1,inod)
114 yi=x(2,inod)
115 zi=x(3,inod)
116 vxi=v(1,inod)
117 vyi=v(2,inod)
118 vzi=v(3,inod)
119 di =spbuf(1,n)
120 rhoi =spbuf(2,n)
121C
122 sxx=wa(1,n)
123 syy=wa(2,n)
124 szz=wa(3,n)
125 sxy=wa(4,n)
126 syz=wa(5,n)
127 sxz=wa(6,n)
128 iprt =ipartsp(n)
129 iprop=ipart(2,iprt)
130 qa = geo(14,iprop)
131 qb = geo(15,iprop)
132C
133C for artificial viscosity computation.
134 divvi=abs(wa(13,n))
135 rotvi=wa(14,n)
136 fi =divvi/max(em20,divvi+rotvi)
137C
138 alphai=wacomp(1,n)
139C BETAXI=WACOMP(2,N)
140C BETAYI=WACOMP(3,N)
141C BETAZI=WACOMP(4,N)
142 alphaxi=wacomp( 5,n)
143 alphayi=wacomp( 6,n)
144 alphazi=wacomp( 7,n)
145 betaxxi=wacomp( 8,n)
146 betayxi=wacomp( 9,n)
147 betazxi=wacomp(10,n)
148 betaxyi=wacomp(11,n)
149 betayyi=wacomp(12,n)
150 betazyi=wacomp(13,n)
151 betaxzi=wacomp(14,n)
152 betayzi=wacomp(15,n)
153 betazzi=wacomp(16,n)
154C----------------------------------
155C FORCES + ARTIFICIAL VISCOUS FORCES + STABILITY:
156C----------------------------------
157 nvois=kxsp(4,n)
158 DO j=1,nvois
159 jnod=ixsp(j,n)
160 IF(jnod>0)THEN
161 m=nod2sp(jnod)
162C
163C Solids to SPH, no interaction if both particles are inactive
164 IF(kxsp(2,n)<=0.AND.kxsp(2,m)<=0)cycle
165C
166 xj=x(1,jnod)
167 yj=x(2,jnod)
168 zj=x(3,jnod)
169 vxj=v(1,jnod)
170 vyj=v(2,jnod)
171 vzj=v(3,jnod)
172 dj =spbuf(1,m)
173 dij=half*(di+dj)
174 rhoj =spbuf(2,m)
175 txx=wa(1,m)
176 tyy=wa(2,m)
177 tzz=wa(3,m)
178 txy=wa(4,m)
179 tyz=wa(5,m)
180 txz=wa(6,m)
181 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
182 wgrdx=wgrad(1)
183 wgrdy=wgrad(2)
184 wgrdz=wgrad(3)
185 wgrad(1)=wgrdx*alphai+wght*alphaxi
186 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
187 wgrad(2)=wgrdy*alphai+wght*alphayi
188 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
189 wgrad(3)=wgrdz*alphai+wght*alphazi
190 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
191! Old order1 correction
192! BETAX=1.+BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
193! WGRAD(1)=WGRDX*ALPHAI*BETAX
194! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
195! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
196! WGRAD(2)=WGRDY*ALPHAI*BETAX
197! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
198! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
199! wgrad(3)=wgrdz*alphai*betax
200! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
201! . BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
202C----------
203C noyau conjugue Grad[Wa(b)]
204 alphaj=wacomp(1,m)
205C BETAXJ=WACOMP(2,M)
206C BETAYJ=WACOMP(3,M)
207C BETAZJ=WACOMP(4,M)
208 alphaxj=wacomp( 5,m)
209 alphayj=wacomp( 6,m)
210 alphazj=wacomp( 7,m)
211 betaxxj=wacomp( 8,m)
212 betayxj=wacomp( 9,m)
213 betazxj=wacomp(10,m)
214 betaxyj=wacomp(11,m)
215 betayyj=wacomp(12,m)
216 betazyj=wacomp(13,m)
217 betaxzj=wacomp(14,m)
218 betayzj=wacomp(15,m)
219 betazzj=wacomp(16,m)
220C
221 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
222 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
223 wgrd(2)=-wgrdy*alphaj+wght*alphayj
224 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
225 wgrd(3)=-wgrdz*alphaj+wght*alphazj
226 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
227C
228! Old order1 correction
229! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
230! WGRD(1)=-WGRDX*ALPHAJ*BETAX
231! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
232! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
233! WGRD(2)=-WGRDY*ALPHAJ*BETAX
234! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
235! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
236! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
237! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
238! . BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
239C
240 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
241 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
242 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
243C BX=TXX*WGRAD(1)+TXY*WGRAD(2)+TXZ*WGRAD(3)
244C BY=TXY*WGRAD(1)+TYY*WGRAD(2)+TYZ*WGRAD(3)
245C BZ=TXZ*WGRAD(1)+TYZ*WGRAD(2)+TZZ*WGRAD(3)
246 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
247 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
248 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
249 mm=spbuf(12,n)*spbuf(12,m)
250C--------
251 vi =spbuf(12,n)/max(em20,rhoi)
252 vj =spbuf(12,m)/max(em20,rhoj)
253 vij=vi*vj
254 fx=vij*(ax+bx)
255 fy=vij*(ay+by)
256 fz=vij*(az+bz)
257C--------
258 wi=zero
259 IF(stab(7,n)/=zero.AND.stab(7,m)/=zero)THEN
260 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
261 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
262 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
263 dx=-(stab(1,m)*wgrd(1)+stab(4,m)*wgrd(2)+stab(6,m)*wgrd(3))
264 dy=-(stab(4,m)*wgrd(1)+stab(2,m)*wgrd(2)+stab(5,m)*wgrd(3))
265 dz=-(stab(6,m)*wgrd(1)+stab(5,m)*wgrd(2)+stab(3,m)*wgrd(3))
266C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
267 ww=wght*dij*dij*dij
268 wr=half*(stab(7,n)+stab(7,m))
269 wi=ww*ww*ww*ww*wr
270 tx=vij*wi*(cx+dx)
271 ty=vij*wi*(cy+dy)
272 tz=vij*wi*(cz+dz)
273 fx=fx+tx
274 fy=fy+ty
275 fz=fz+tz
276 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
277 END IF
278C--------
279C artitifial viscosity.
280 dxij=(vxi-vxj)*(xi-xj)
281 . +(vyi-vyj)*(yi-yj)
282 . +(vzi-vzj)*(zi-zj)
283 nxij=(xi-xj)*(xi-xj)
284 . +(yi-yj)*(yi-yj)
285 . +(zi-zj)*(zi-zj)
286 muij=dij*dxij/(nxij+em02*dij*dij)
287 divvj=abs(wa(13,m))
288 rotvj=wa(14,m)
289 fj =divvj/max(em20,divvj+rotvj)
290 muij=min(muij,zero)
291 muij=muij*(fi+fj)*half
292 ssp=(wa(8,n)+wa(8,m))*half
293 muij2=muij*muij
294 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
295 fact=mm*pij
296C---------
297C FV(j,i)=-FV(i,j)
298 wgrdx=(wgrad(1)-wgrd(1))*half
299 wgrdy=(wgrad(2)-wgrd(2))*half
300 wgrdz=(wgrad(3)-wgrd(3))*half
301 fvx=-fact*wgrdx
302 fvy=-fact*wgrdy
303 fvz=-fact*wgrdz
304C
305 IF((nodadt/=0).OR.(i7kglo/=0))THEN
306C
307C nodal stability
308 dldt=abs(dxij)
309 l=sqrt(nxij)
310 dldt=max(em20,dldt/max(em20,l))
311C
312 volj=spbuf(12,m)/max(em20,rhoj)
313 drhoidr= volj
314 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
315 ssp2i=wa(9,n)*wa(9,n)
316 stii = volj*spbuf(12,n)*ssp2i*drhoidr
317C
318 voli =spbuf(12,n)/max(em20,rhoi)
319 drhojdr= voli
320 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
321 ssp2j=wa(9,m)*wa(9,m)
322 stij = voli*spbuf(12,m)*ssp2j*drhojdr
323C
324 stij=two*(stii+stij)
325C
326C K*= fv divided by dl (c = fv divided by v)
327 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
328 cij =mm*abs(pij)*wnorm/dldt
329 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
330 sfac =sqrt(1.+dzeta*dzeta)-dzeta
331 stij =stij/max(em20,sfac*sfac)
332 wa(7,n)=wa(7,n)+stij*(one+wi)
333 ENDIF
334 ELSE ! cellule remote
335 nn = -jnod
336C
337C Solids to SPH, no interaction if both particles are inactive
338 IF(kxsp(2,n)<=0.AND.xsphr(13,nn)<=0)cycle
339 xj=xsphr(3,nn)
340 yj=xsphr(4,nn)
341 zj=xsphr(5,nn)
342 vxj=xsphr(9,nn)
343 vyj=xsphr(10,nn)
344 vzj=xsphr(11,nn)
345 dj =xsphr(2,nn)
346 dij=half*(di+dj)
347 rhoj =xsphr(7,nn)
348 txx=war(1,nn)
349 tyy=war(2,nn)
350 tzz=war(3,nn)
351 txy=war(4,nn)
352 tyz=war(5,nn)
353 txz=war(6,nn)
354 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
355 wgrdx=wgrad(1)
356 wgrdy=wgrad(2)
357 wgrdz=wgrad(3)
358 wgrad(1)=wgrdx*alphai+wght*alphaxi
359 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
360 wgrad(2)=wgrdy*alphai+wght*alphayi
361 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
362 wgrad(3)=wgrdz*alphai+wght*alphazi
363 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
364! Old order1 correction
365! BETAX=1.+BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
366! WGRAD(1)=WGRDX*ALPHAI*BETAX
367! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
368! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
369! WGRAD(2)=WGRDY*ALPHAI*BETAX
370! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
371! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
372! WGRAD(3)=WGRDZ*ALPHAI*BETAX
373! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
374! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
375C----------
376C noyau conjugue Grad[Wa(b)]
377 alphaj=wacompr(1,nn)
378C BETAXJ=WACOMPR(2,NN)
379C BETAYJ=WACOMPR(3,NN)
380C BETAZJ=WACOMPR(4,NN)
381 alphaxj=wacompr( 5,nn)
382 alphayj=wacompr( 6,nn)
383 alphazj=wacompr( 7,nn)
384 betaxxj=wacompr( 8,nn)
385 betayxj=wacompr( 9,nn)
386 betazxj=wacompr(10,nn)
387 betaxyj=wacompr(11,nn)
388 betayyj=wacompr(12,nn)
389 betazyj=wacompr(13,nn)
390 betaxzj=wacompr(14,nn)
391 betayzj=wacompr(15,nn)
392 betazzj=wacompr(16,nn)
393C
394 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
395 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
396 wgrd(2)=-wgrdy*alphaj+wght*alphayj
397 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
398 wgrd(3)=-wgrdz*alphaj+wght*alphazj
399 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
400! Old order1 correction
401! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
402! WGRD(1)=-WGRDX*ALPHAJ*BETAX
403! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
404! . (BETAXXJ*(XJ-XI)+BETAYXJ*(YJ-YI)+BETAZXJ*(ZJ-ZI)+BETAXJ))
405! WGRD(2)=-WGRDY*ALPHAJ*BETAX
406! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
407! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
408! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
409! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
410! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
411C
412 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
413 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
414 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
415C BX=TXX*WGRAD(1)+TXY*WGRAD(2)+TXZ*WGRAD(3)
416C BY=TXY*WGRAD(1)+TYY*WGRAD(2)+TYZ*WGRAD(3)
417C BZ=TXZ*WGRAD(1)+TYZ*WGRAD(2)+TZZ*WGRAD(3)
418 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
419 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
420 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
421 mm=spbuf(12,n)*xsphr(8,nn)
422C--------
423 vi =spbuf(12,n)/max(em20,rhoi)
424 vj =xsphr(8,nn)/max(em20,rhoj)
425 vij=vi*vj
426 fx=vij*(ax+bx)
427 fy=vij*(ay+by)
428 fz=vij*(az+bz)
429C--------
430 wi=zero
431 IF(stab(7,n)/=zero.AND.stab(7,numsph+nn)/=zero)THEN
432 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
433 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
434 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
435 nr=numsph+nn
436 dx=-(stab(1,nr)*wgrd(1)+stab(4,nr)*wgrd(2)+stab(6,nr)*wgrd(3))
437 dy=-(stab(4,nr)*wgrd(1)+stab(2,nr)*wgrd(2)+stab(5,nr)*wgrd(3))
438 dz=-(stab(6,nr)*wgrd(1)+stab(5,nr)*wgrd(2)+stab(3,nr)*wgrd(3))
439C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
440 ww=wght*dij*dij*dij
441 wr=half*(stab(7,n)+stab(7,numsph+nn))
442 wi=ww*ww*ww*ww*wr
443 tx=vij*wi*(cx+dx)
444 ty=vij*wi*(cy+dy)
445 tz=vij*wi*(cz+dz)
446 fx=fx+tx
447 fy=fy+ty
448 fz=fz+tz
449 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
450 END IF
451C--------
452C artitifial viscosity.
453 dxij=(vxi-vxj)*(xi-xj)
454 . +(vyi-vyj)*(yi-yj)
455 . +(vzi-vzj)*(zi-zj)
456 nxij=(xi-xj)*(xi-xj)
457 . +(yi-yj)*(yi-yj)
458 . +(zi-zj)*(zi-zj)
459 muij=dij*dxij/(nxij+em02*dij*dij)
460C WA(13) <=> WAR(9) ; WA(14) <=> WAR(10)
461 divvj=abs(war(9,nn))
462 rotvj=war(10,nn)
463 fj =divvj/max(em20,divvj+rotvj)
464 muij=min(muij,zero)
465 muij=muij*(fi+fj)*half
466C WA(8) <=> WAR(7)
467 ssp=(wa(8,n)+war(7,nn))*half
468 muij2=muij*muij
469 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
470 fact=mm*pij
471C---------
472C FV(j,i)=-FV(i,j)
473 wgrdx=(wgrad(1)-wgrd(1))*half
474 wgrdy=(wgrad(2)-wgrd(2))*half
475 wgrdz=(wgrad(3)-wgrd(3))*half
476 fvx=-fact*wgrdx
477 fvy=-fact*wgrdy
478 fvz=-fact*wgrdz
479 IF((nodadt/=0).OR.(i7kglo/=0))THEN
480C
481C nodal stability
482 dldt=abs(dxij)
483 l=sqrt(nxij)
484 dldt=max(em20,dldt/max(em20,l))
485C
486 volj=xsphr(8,nn)/max(em20,rhoj)
487 drhoidr= volj
488 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
489 ssp2i=wa(9,n)*wa(9,n)
490 stii = volj*spbuf(12,n)*ssp2i*drhoidr
491C
492 voli =spbuf(12,n)/max(em20,rhoi)
493 drhojdr= voli
494 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
495C WA(9) <=> WAR(8)
496 ssp2j=war(8,nn)*war(8,nn)
497 stij = voli*xsphr(8,nn)*ssp2j*drhojdr
498C
499 stij=two*(stii+stij)
500C
501C K*= fv divided by dl (c = fv divided by v)
502 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
503 cij =mm*abs(pij)*wnorm/dldt
504 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
505 sfac =sqrt(1.+dzeta*dzeta)-dzeta
506 stij =stij/max(em20,sfac*sfac)
507 wa(7,n)=wa(7,n)+stij*(one+wi)
508 ENDIF
509 END IF
510C
511 fx=fx+fvx
512 fy=fy+fvy
513 fz=fz+fvz
514 wa(10,n)=wa(10,n)+fx
515 wa(11,n)=wa(11,n)+fy
516 wa(12,n)=wa(12,n)+fz
517C
518C for work of artificial viscosity forces (per cell).
519 wvis = min(zero,
520 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
521 spbuf(11,n)=spbuf(11,n)+wvis
522 ENDDO
523C----------------------------------
524 nvoiss=kxsp(6,n)
525 DO 200 j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
526 js=ixsp(j,n)
527 IF(js>0)THEN
528 sm=js/(nspcond+1)
529C
530C Solids to SPH, no interaction if both particles are inactive
531 IF(kxsp(2,n)<=0.AND.kxsp(2,sm)<=0)cycle
532C
533 nc=mod(js,nspcond+1)
534 js=ispsym(nc,sm)
535 xj =xspsym(1,js)
536 yj =xspsym(2,js)
537 zj =xspsym(3,js)
538 vxj=vspsym(1,js)
539 vyj=vspsym(2,js)
540 vzj=vspsym(3,js)
541 dj =spbuf(1,sm)
542 dij =half*(di+dj)
543 rhoj=spbuf(2,sm)
544 jnod=kxsp(3,sm)
545C
546 txx=wasigsm(1,js)
547 tyy=wasigsm(2,js)
548 tzz=wasigsm(3,js)
549 txy=wasigsm(4,js)
550 tyz=wasigsm(5,js)
551 txz=wasigsm(6,js)
552 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
553 wgrdx=wgrad(1)
554 wgrdy=wgrad(2)
555 wgrdz=wgrad(3)
556 wgrad(1)=wgrdx*alphai+wght*alphaxi
557 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
558 wgrad(2)=wgrdy*alphai+wght*alphayi
559 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
560 wgrad(3)=wgrdz*alphai+wght*alphazi
561 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
562! Old order1 correction
563! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
564! WGRAD(1)=WGRDX*ALPHAI*BETAX
565! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
566! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
567! WGRAD(2)=WGRDY*ALPHAI*BETAX
568! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
569! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
570! WGRAD(3)=WGRDZ*ALPHAI*BETAX
571! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
572! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
573C----------
574C nucleus combines.
575 alphaj=wacomp(1,sm)
576C BETAXJ=WACOMP(2,SM)
577C BETAYJ=WACOMP(3,SM)
578C BETAZJ=WACOMP(4,SM)
579 betaxj=wsmcomp(1,js)
580 betayj=wsmcomp(2,js)
581 betazj=wsmcomp(3,js)
582C ALPHAXJ=WACOMP( 5,SM)
583C ALPHAYJ=WACOMP( 6,SM)
584C ALPHAZJ=WACOMP( 7,SM)
585 alphaxj=wsmcomp( 4,js)
586 alphayj=wsmcomp( 5,js)
587 alphazj=wsmcomp( 6,js)
588 betaxxj=wacomp( 8,sm)
589 betayxj=wacomp( 9,sm)
590 betazxj=wacomp(10,sm)
591 betaxyj=wacomp(11,sm)
592 betayyj=wacomp(12,sm)
593 betazyj=wacomp(13,sm)
594 betaxzj=wacomp(14,sm)
595 betayzj=wacomp(15,sm)
596 betazzj=wacomp(16,sm)
597 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
598 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
599 wgrd(2)=-wgrdy*alphaj+wght*alphayj
600 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
601 wgrd(3)=-wgrdz*alphaj+wght*alphazj
602 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
603! Old order1 correction
604! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
605! WGRD(1)=-WGRDX*ALPHAJ*BETAX
606! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
607! . (betaxxj*(xj-xi)+betayxj*(yj-yi)+betazxj*(zj-zi)+betaxj))
608! wgrd(2)=-wgrdy*alphaj*betax
609! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
610! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
611! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
612! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
613! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
614 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
615 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
616 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
617 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
618 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
619 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
620 mm=spbuf(12,n)*spbuf(12,sm)
621C--------
622 vi =spbuf(12,n)/max(em20,rhoi)
623 vj =spbuf(12,sm)/max(em20,rhoj)
624 vij=vi*vj
625 fx=vij*(ax+bx)
626 fy=vij*(ay+by)
627 fz=vij*(az+bz)
628 wi=zero
629 IF(stab(7,n)/=zero.AND.stab(7,sm)/=zero)THEN
630 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
631 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
632 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
633 ks=numsph+nsphr+js
634 dx=-(stab(1,ks)*wgrd(1)+stab(4,ks)*wgrd(2)+stab(6,ks)*wgrd(3))
635 dy=-(stab(4,ks)*wgrd(1)+stab(2,ks)*wgrd(2)+stab(5,ks)*wgrd(3))
636 dz=-(stab(6,ks)*wgrd(1)+stab(5,ks)*wgrd(2)+stab(3,ks)*wgrd(3))
637C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
638 ww=wght*dij*dij*dij
639 wr=half*(stab(7,n)+stab(7,sm))
640 wi=ww*ww*ww*ww*wr
641 tx=vij*wi*(cx+dx)
642 ty=vij*wi*(cy+dy)
643 tz=vij*wi*(cz+dz)
644 fx=fx+tx
645 fy=fy+ty
646 fz=fz+tz
647 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
648 END IF
649C--------
650C artitifial viscosity.
651 dxij=(vxi-vxj)*(xi-xj)
652 . +(vyi-vyj)*(yi-yj)
653 . +(vzi-vzj)*(zi-zj)
654 nxij=(xi-xj)*(xi-xj)
655 . +(yi-yj)*(yi-yj)
656 . +(zi-zj)*(zi-zj)
657 muij=dij*dxij/(nxij+em02*dij*dij)
658 divvj=abs(wa(13,sm))
659 rotvj=wa(14,sm)
660 fj=divvj/max(em20,divvj+rotvj)
661 muij=min(muij,zero)
662 muij=muij*(fi+fj)*half
663 ssp=(wa(8,n)+wa(8,sm))*half
664 muij2=muij*muij
665 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
666 fact=mm*pij
667C---------
668 wgrdx=(wgrad(1)-wgrd(1))*half
669 wgrdy=(wgrad(2)-wgrd(2))*half
670 wgrdz=(wgrad(3)-wgrd(3))*half
671 fvx=-fact*wgrdx
672 fvy=-fact*wgrdy
673 fvz=-fact*wgrdz
674 IF((nodadt/=0).OR.(i7kglo/=0))THEN
675C nodal stability
676 dldt=abs(dxij)
677 l =sqrt(nxij)
678 dldt=max(em20,dldt/max(em20,l))
679C
680 volj=spbuf(12,sm)/max(em20,rhoj)
681 drhoidr= volj
682 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
683 ssp2i=wa(9,n)*wa(9,n)
684 stii = volj*spbuf(12,n)*ssp2i*drhoidr
685C
686 voli =spbuf(12,n)/max(em20,rhoi)
687 drhojdr= voli
688 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
689 ssp2j=wa(9,sm)*wa(9,sm)
690 stij = voli*spbuf(12,sm)*ssp2j*drhojdr
691C
692 stij=two*(stii+stij)
693C
694C K*= fv divided by dl (c = fv divided by v)
695 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
696 cij =mm*abs(pij)*wnorm/dldt
697 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
698 sfac =sqrt(one +dzeta*dzeta)-dzeta
699 stij =stij/max(em20,sfac*sfac)
700 wa(7,n)=wa(7,n)+stij*(one+wi)
701 ENDIF
702 fx=fx+fvx
703 fy=fy+fvy
704 fz=fz+fvz
705 wa(10,n)=wa(10,n)+fx
706 wa(11,n)=wa(11,n)+fy
707 wa(12,n)=wa(12,n)+fz
708C
709C for work of artificial viscosity forces (per cell).
710 wvis = min(zero,
711 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
712 spbuf(11,n)=spbuf(11,n)+wvis
713 ELSE ! Symmetrical particle particle Remote
714 sm=-js/(nspcond+1)
715C
716C Solids to SPH, no interaction if both particles are inactive
717 IF(kxsp(2,n)<=0.AND.xsphr(13,sm)<=0)cycle
718 nc=mod(-js,nspcond+1)
719 js=ispsymr(nc,sm)
720 xj =xspsym(1,js)
721 yj =xspsym(2,js)
722 zj =xspsym(3,js)
723 vxj=vspsym(1,js)
724 vyj=vspsym(2,js)
725 vzj=vspsym(3,js)
726 dj =xsphr(2,sm)
727 dij =half*(di+dj)
728 rhoj=xsphr(7,sm)
729C
730 txx=wasigsm(1,js)
731 tyy=wasigsm(2,js)
732 tzz=wasigsm(3,js)
733 txy=wasigsm(4,js)
734 tyz=wasigsm(5,js)
735 txz=wasigsm(6,js)
736 CALL weight1(xi,yi,zi,xj,yj,zj,dij,wght,wgrad)
737 wgrdx=wgrad(1)
738 wgrdy=wgrad(2)
739 wgrdz=wgrad(3)
740 wgrad(1)=wgrdx*alphai+wght*alphaxi
741 . +wgrdx*betaxxi+wgrdy*betaxyi+wgrdz*betaxzi
742 wgrad(2)=wgrdy*alphai+wght*alphayi
743 . +wgrdx*betayxi+wgrdy*betayyi+wgrdz*betayzi
744 wgrad(3)=wgrdz*alphai+wght*alphazi
745 . +wgrdx*betazxi+wgrdy*betazyi+wgrdz*betazzi
746! Old order1 correction
747! BETAX=ONE +BETAXI*(XI-XJ)+BETAYI*(YI-YJ)+BETAZI*(ZI-ZJ)
748! WGRAD(1)=WGRDX*ALPHAI*BETAX
749! . +WGHT*(ALPHAXI*BETAX+ALPHAI*
750! . (BETAXXI*(XI-XJ)+BETAYXI*(YI-YJ)+BETAZXI*(ZI-ZJ)+BETAXI))
751! WGRAD(2)=WGRDY*ALPHAI*BETAX
752! . +WGHT*(ALPHAYI*BETAX+ALPHAI*
753! . (BETAXYI*(XI-XJ)+BETAYYI*(YI-YJ)+BETAZYI*(ZI-ZJ)+BETAYI))
754! WGRAD(3)=WGRDZ*ALPHAI*BETAX
755! . +WGHT*(ALPHAZI*BETAX+ALPHAI*
756! . (BETAXZI*(XI-XJ)+BETAYZI*(YI-YJ)+BETAZZI*(ZI-ZJ)+BETAZI))
757C----------
758C nucleus combines.
759 alphaj=wacompr(1,sm)
760C BETAXJ=WACOMPR(2,SM)
761C BETAYJ=WACOMPR(3,SM)
762C BETAZJ=WACOMPR(4,SM)
763 betaxj=wsmcomp(1,js)
764 betayj=wsmcomp(2,js)
765 betazj=wsmcomp(3,js)
766C ALPHAXJ=WACOMPR( 5,SM)
767C ALPHAYJ=WACOMPR( 6,SM)
768C ALPHAZJ=WACOMPR( 7,SM)
769 alphaxj=wsmcomp( 4,js)
770 alphayj=wsmcomp( 5,js)
771 alphazj=wsmcomp( 6,js)
772 betaxxj=wacompr( 8,sm)
773 betayxj=wacompr( 9,sm)
774 betazxj=wacompr(10,sm)
775 betaxyj=wacompr(11,sm)
776 betayyj=wacompr(12,sm)
777 betazyj=wacompr(13,sm)
778 betaxzj=wacompr(14,sm)
779 betayzj=wacompr(15,sm)
780 betazzj=wacompr(16,sm)
781C
782 wgrd(1)=-wgrdx*alphaj+wght*alphaxj
783 . -wgrdx*betaxxj-wgrdy*betaxyj-wgrdz*betaxzj
784 wgrd(2)=-wgrdy*alphaj+wght*alphayj
785 . -wgrdx*betayxj-wgrdy*betayyj-wgrdz*betayzj
786 wgrd(3)=-wgrdz*alphaj+wght*alphazj
787 . -wgrdx*betazxj-wgrdy*betazyj-wgrdz*betazzj
788! Old order1 correction
789! BETAX=ONE +BETAXJ*(XJ-XI)+BETAYJ*(YJ-YI)+BETAZJ*(ZJ-ZI)
790! WGRD(1)=-WGRDX*ALPHAJ*BETAX
791! . +WGHT*(ALPHAXJ*BETAX+ALPHAJ*
792! . (betaxxj*(xj-xi)+betayxj*(yj-yi)+betazxj*(zj-zi)+betaxj))
793! WGRD(2)=-WGRDY*ALPHAJ*BETAX
794! . +WGHT*(ALPHAYJ*BETAX+ALPHAJ*
795! . (BETAXYJ*(XJ-XI)+BETAYYJ*(YJ-YI)+BETAZYJ*(ZJ-ZI)+BETAYJ))
796! WGRD(3)=-WGRDZ*ALPHAJ*BETAX
797! . +WGHT*(ALPHAZJ*BETAX+ALPHAJ*
798! . (BETAXZJ*(XJ-XI)+BETAYZJ*(YJ-YI)+BETAZZJ*(ZJ-ZI)+BETAZJ))
799 ax=sxx*wgrad(1)+sxy*wgrad(2)+sxz*wgrad(3)
800 ay=sxy*wgrad(1)+syy*wgrad(2)+syz*wgrad(3)
801 az=sxz*wgrad(1)+syz*wgrad(2)+szz*wgrad(3)
802 bx=-(txx*wgrd(1)+txy*wgrd(2)+txz*wgrd(3))
803 by=-(txy*wgrd(1)+tyy*wgrd(2)+tyz*wgrd(3))
804 bz=-(txz*wgrd(1)+tyz*wgrd(2)+tzz*wgrd(3))
805 mm=spbuf(12,n)*xsphr(8,sm)
806C--------
807 vi =spbuf(12,n)/max(em20,rhoi)
808 vj =xsphr(8,sm)/max(em20,rhoj)
809 vij=vi*vj
810 fx=vij*(ax+bx)
811 fy=vij*(ay+by)
812 fz=vij*(az+bz)
813 wi=zero
814 IF(stab(7,n)/=zero.AND.stab(7,numsph+sm)/=zero)THEN
815 cx=stab(1,n)*wgrad(1)+stab(4,n)*wgrad(2)+stab(6,n)*wgrad(3)
816 cy=stab(4,n)*wgrad(1)+stab(2,n)*wgrad(2)+stab(5,n)*wgrad(3)
817 cz=stab(6,n)*wgrad(1)+stab(5,n)*wgrad(2)+stab(3,n)*wgrad(3)
818 ks=numsph+nsphr+js
819 dx=-(stab(1,ks)*wgrd(1)+stab(4,ks)*wgrd(2)+stab(6,ks)*wgrd(3))
820 dy=-(stab(4,ks)*wgrd(1)+stab(2,ks)*wgrd(2)+stab(5,ks)*wgrd(3))
821 dz=-(stab(6,ks)*wgrd(1)+stab(5,ks)*wgrd(2)+stab(3,ks)*wgrd(3))
822C CALL WEIGHT0(XI,YI,ZI,XJ,YJ,ZJ,DIJ,WGHT)
823 ww=wght*dij*dij*dij
824 wr=half*(stab(7,n)+stab(7,numsph+sm))
825 wi=ww*ww*ww*ww*wr
826 tx=vij*wi*(cx+dx)
827 ty=vij*wi*(cy+dy)
828 tz=vij*wi*(cz+dz)
829 fx=fx+tx
830 fy=fy+ty
831 fz=fz+tz
832 wfextt=wfextt+(tx*v(1,inod)+ty*v(2,inod)+tz*v(3,inod))*dt1
833 END IF
834C--------
835C artitifial viscosity.
836 dxij=(vxi-vxj)*(xi-xj)
837 . +(vyi-vyj)*(yi-yj)
838 . +(vzi-vzj)*(zi-zj)
839 nxij=(xi-xj)*(xi-xj)
840 . +(yi-yj)*(yi-yj)
841 . +(zi-zj)*(zi-zj)
842 muij=dij*dxij/(nxij+em02*dij*dij)
843C WA(13) <=> WAR(9) ; WA(14) <=> WAR(10)
844 divvj=abs(war(9,sm))
845 rotvj=war(10,sm)
846 fj=divvj/max(em20,divvj+rotvj)
847 muij=min(muij,zero)
848 muij=muij*(fi+fj)*half
849C WA(8) <=> WAR(7)
850 ssp=(wa(8,n)+war(7,sm))*half
851 muij2=muij*muij
852 pij =(qa*muij2-qb*ssp*muij)*two/max(em20,rhoi+rhoj)
853 fact=mm*pij
854C---------
855 wgrdx=(wgrad(1)-wgrd(1))*half
856 wgrdy=(wgrad(2)-wgrd(2))*half
857 wgrdz=(wgrad(3)-wgrd(3))*half
858 fvx=-fact*wgrdx
859 fvy=-fact*wgrdy
860 fvz=-fact*wgrdz
861 IF((nodadt/=0).OR.(i7kglo/=0))THEN
862C nodal stability
863 dldt=abs(dxij)
864 l =sqrt(nxij)
865 dldt=max(em20,dldt/max(em20,l))
866C
867 volj=xsphr(8,sm)/max(em20,rhoj)
868 drhoidr= volj
869 . * (wgrad(1)*wgrad(1)+wgrad(2)*wgrad(2)+wgrad(3)*wgrad(3))
870 ssp2i=wa(9,n)*wa(9,n)
871 stii = volj*spbuf(12,n)*ssp2i*drhoidr
872C
873 voli =spbuf(12,n)/max(em20,rhoi)
874 drhojdr= voli
875 . * (wgrd(1)*wgrd(1)+wgrd(2)*wgrd(2)+wgrd(3)*wgrd(3))
876C WA(9) <=> WAR(8)
877 ssp2j=war(8,sm)*war(8,sm)
878 stij = voli*xsphr(8,sm)*ssp2j*drhojdr
879C
880 stij=two*(stii+stij)
881C
882C K*= fv divided by dl (c = fv divided by v)
883 wnorm=sqrt(wgrdx*wgrdx+wgrdy*wgrdy+wgrdz*wgrdz)
884 cij =mm*abs(pij)*wnorm/dldt
885 dzeta=cij/max(em20,sqrt(two*stij*spbuf(12,n)))
886 sfac =sqrt(one +dzeta*dzeta)-dzeta
887 stij =stij/max(em20,sfac*sfac)
888 wa(7,n)=wa(7,n)+stij*(one+wi)
889 ENDIF
890 fx=fx+fvx
891 fy=fy+fvy
892 fz=fz+fvz
893 wa(10,n)=wa(10,n)+fx
894 wa(11,n)=wa(11,n)+fy
895 wa(12,n)=wa(12,n)+fz
896C
897C for work of artificial viscosity forces (per cell).
898 wvis = min(zero,
899 . half*(fvx*(vxi-vxj)+fvy*(vyi-vyj)+fvz*(vzi-vzj)))
900 spbuf(11,n)=spbuf(11,n)+wvis
901 END IF
902 200 CONTINUE
903C-------
904 10 CONTINUE
905#include "lockon.inc"
906 wfext=wfext+wfextt
907#include "lockoff.inc"
908C----------------------------------
909 RETURN
910 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
spforcp
subroutine spforcp(pm, geo, x, v, ms, spbuf, itab, pld, bufmat, bufgeo, partsav, fsav, dt2t, iparg, npc, kxsp, ixsp, nod2sp, neltst, ityptst, ipart, ipartsp, ispcond, xframe, ispsym, xspsym, vspsym, wa, wasigsm, wacomp, wsmcomp, waspact, war, stab, wfext)
Definition spforcp.F:40
weight1
subroutine weight1(xi, yi, zi, xj, yj, zj, h, w, wgrad)
Definition weight.F:79