source: CIVL/examples/fortran/nek5000/core/cmt/diagnostics.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 19.2 KB
Line 
1 subroutine printalot(e,f,lu,pileoffaces)
2 include 'SIZE'
3 include 'DG'
4 include 'CMTDATA'
5 include 'GEOM'
6 parameter (ldd=lxd*lyd*lzd)
7 common /ctmp1/ xface(lx1*lz1,2*ldim),yface(lx1*lz1,2*ldim),
8 > supapad(6*ldd-4*ldim*lx1*lz1)
9 real pileoffaces(lx1*lz1,2*ldim,nelt,toteq)!,ldim)
10 integer e,f,lu
11 real ghere, betahere,pi,kond,h
12 pi=4.0*atan(1.0)
13 ghere=1.4
14 betahere=5.0
15
16 call full2face_cmt(1,lx1,ly1,lz1,iface_flux(1,e),xface,
17 > xm1(1,1,1,e))
18 call full2face_cmt(1,lx1,ly1,lz1,iface_flux(1,e),yface,
19 > ym1(1,1,1,e))
20 nxz=lx1*lz1
21 write(lu,*) '# e=',e
22 do i=1,nxz
23! x=xface(i,f)-5.0
24! y=yface(i,f)
25! r2=x**2+y**2
26! zeu=-betahere*y*0.5/pi*exp(1.0-r2)
27! zev= betahere*x*0.5/pi*exp(1.0-r2)
28! tau11=2.0*muref*betahere*x*y/pi*exp(1.0-r2)
29! tau12=betahere/pi*(y**2-x**2)*exp(1.0-r2)*muref
30! tau22=-tau11
31! zework1=zeu*tau11+zev*tau12
32! zework2=zeu*tau12+zev*tau22
33! kond=cpgref*muref/prlam
34! tx=0.25*betahere**2/pi/pi*(ghere-1.0)/ghere*exp(2.*(1.0-r2))*x
35! ty=0.25*betahere**2/pi/pi*(ghere-1.0)/ghere*exp(2.*(1.0-r2))*y
36! fluxx=zework1-kond*tx
37! fluxy=zework2-kond*ty
38 y=yface(i,f)
39 u2=1.0
40 h=1.0
41 tau12=u2*muref/h
42! write(lu,'(5e17.5)') xface(i,f),yface(i,f),-tau12,
43! > pileoffaces(i,f,e,3,1),
44! > pileoffaces(i,f,e,3,2)
45! > pileoffaces(i,f,e,3,1),-tau11,
46! write(lu,'(4e17.5)') xface(i,f),yface(i,f),
47! > pileoffaces(i,f,e,3,1),
48! > pileoffaces(i,f,e,3,2)
49 write(lu,*) xface(i,f),yface(i,f), pileoffaces(i,f,e,1)
50 enddo
51
52 return
53 end
54c-----------------------------------------------------------------------
55! A bunch of diagostic routines
56c-----------------------------------------------------------------------
57 subroutine matout_rowsum(ab,na,nb)
58 integer na,nb
59 real ab(na,nb)
60 character*100 zefmt
61 write(zefmt,'(a1,i2,a6)') '(',nb,'e15.7)'
62 do i=1,na
63 write(6,zefmt) (ab(i,j),j=1,nb)
64! write(6,*) i,'rowsum=',sum(ab(i,1:nb))
65 enddo
66 return
67 end
68c----------------------------------------------------------------------
69 subroutine out_fld_nek
70 include 'SIZE'
71 include 'SOLN'
72 COMMON /solnconsvar/ U(LX1,LY1,LZ1,TOTEQ,lelt)
73 COMMON /SCRNS/ OTVAR(LX1,LY1,LZ1,lelt,7)
74 real OTVAR
75 real phig(lx1,ly1,lz1,lelt)
76 common /otherpvar/ phig
77 integer e,f
78
79 n = lx1*ly1*lz1
80 do e=1,nelt
81 call copy(otvar(1,1,1,e,4),u(1,1,1,1,e),n)
82 call copy(otvar(1,1,1,e,5),u(1,1,1,2,e),n)
83 call copy(otvar(1,1,1,e,6),u(1,1,1,3,e),n)
84 call copy(otvar(1,1,1,e,7),u(1,1,1,4,e),n)
85 call copy(otvar(1,1,1,e,1),u(1,1,1,5,e),n)
86 enddo
87
88c call copy(otvar(1,1,1,1,2),tlag(1,1,1,1,1,2),n*nelt) ! s_{n-1}
89c call copy(otvar(1,1,1,1,3),tlag(1,1,1,1,2,1),n*nelt) ! s_n
90 call copy(otvar(1,1,1,1,2),phig(1,1,1,1),n*nelt) ! s_{n-1}
91 call copy(otvar(1,1,1,1,3),pr(1,1,1,1),n*nelt) ! s_n
92
93c ifxyo=.true.
94 if (lx2.ne.lx1) call exitti('Set LX1=LX2 for I/O$',lx2)
95
96 itmp = 3
97 call outpost2(otvar(1,1,1,1,5),otvar(1,1,1,1,6),otvar(1,1,1,1,7)
98 $ ,otvar(1,1,1,1,4),otvar(1,1,1,1,1),itmp,'SLN')
99 return
100 end
101
102c----------------------------------------------------------------------
103
104 subroutine dumpresidue(wfnav,inmbr)
105 include 'SIZE'
106 include 'INPUT'
107 include 'GEOM'
108 include 'MASS'
109 include 'CMTDATA'
110
111 character*32 wfnav
112 character*32 citer,citer2
113 integer e,i1,is,il,i,inmbr,length,eq,is2,il2
114 real rhseqs(toteq)
115
116c write(6,*)wfnav,inmbr
117 i1 =index(wfnav,' ')-1
118 write(citer,*)inmbr
119 write(citer2,*)nid
120 length = len(wfnav)
121 do i=1,length
122 if(citer(i:i).ne.' ')then
123 is=i
124c exit ! Not valid w/ pgf77
125 endif
126 enddo
127 do i=length,1,-1
128 if(citer(i:i).ne.' ')then
129 il=i
130c exit ! Not valid w/ pgf77
131 endif
132 enddo
133! get the string that contains character nid
134 do i=1,length
135 if(citer2(i:i).ne.' ')then
136 is2=i
137c exit ! Not valid w/ pgf77
138 endif
139 enddo
140 do i=length,1,-1
141 if(citer2(i:i).ne.' ')then
142 il2=i
143c exit ! Not valid w/ pgf77
144 endif
145 enddo
146 nxyz1 = lx1*ly1*lz1
147 nxy1 = lx1*ly1
148 open(unit=11,file=wfnav(1:i1)//'.'//'it='//citer(is:il)//
149 $ '.'//'p='//citer2(is2:il2))
150 write(11,*)'Title="',wfnav(1:i1),'"'
151c write(6,*)wfnav(1:i1),'.',citer(is:il)
152 do i=1,length
153 wfnav(i:i)=' '
154 enddo
155 do i=1,length
156 citer(i:i)=' '
157 enddo
158 if(if3d)then
159 write(11,*)'Variables=x,y,z,e1,e2,e3,e4,e5'
160 do e = 1,nelt
161 write(11,*)'zone T="',e,'",i=',lx1,',j=',ly1,',k=',lz1
162 do i=1,nxyz1
163 do eq=1,toteq
164 rhseqs(eq) = res1(i,1,1,e,eq)/bm1(i,1,1,e)
165 enddo
166 write(11,101)xm1(i,1,1,e),ym1(i,1,1,e),zm1(i,1,1,e)
167 $ ,rhseqs(1),rhseqs(2),rhseqs(3),rhseqs(4)
168 $ ,rhseqs(5)
169 enddo
170 enddo
171 else
172 write(11,*)'Variables=x,y,e1,e2,e3,e4,e5'
173 do e = 1,nelt
174 write(11,*)'zone T="',e,'",i=',lx1,',j=',ly1
175 do i=1,nxy1
176 do eq=1,toteq
177 rhseqs(eq) = res1(i,1,1,e,eq)/bm1(i,1,1,e)
178 enddo
179 write(11,102)xm1(i,1,1,e),ym1(i,1,1,e)
180 $ ,rhseqs(1),rhseqs(2),rhseqs(3),rhseqs(4)
181 $ ,rhseqs(5)
182 enddo
183 enddo
184 endif
185 close(11)
186101 format(8(3x,e12.5))
187102 format(7(3x,e14.7))
188 return
189 end
190c----------------------------------------------------------------------
191 subroutine mass_balance(if3d)
192 INCLUDE 'SIZE'
193 INCLUDE 'GEOM'
194 INCLUDE 'MASS'
195 INCLUDE 'TSTEP'
196 COMMON /solnconsvar/ U(LX1,LY1,LZ1,TOTEQ,lelt)
197 logical if3d
198 integer e
199 real msum,total_mass
200c First get the mass in the local domain. Then use
201c Global sum to get the total mass
202c mass = \rho_g \phi_g * Wts(i,j,k)
203c Note - this is only for single processor runs
204c need to add glsum to make it parallel
205 if (if3d)then
206 msum = 0.0
207 do e=1,nelt
208 il = 0
209 do k=1,lz1
210 do j=1,ly1
211 do i=1,lx1
212 il = il +1
213 msum = msum + (u(i,j,k,1,e)*bm1(i,j,k,e))
214 enddo
215 enddo
216 enddo
217 enddo
218 else
219 msum = 0.0
220 do e=1,nelt
221 il = 0
222 do j=1,ly1
223 do i=1,lx1
224 il = il +1
225 msum = msum + (u(i,j,1,1,e)*bm1(i,j,1,e))
226 enddo
227 enddo
228 enddo
229 endif
230 total_mass = glsum(msum,1)
231 if(nio.eq.0)
232c $ write(6,*)'Total mass in the domain ',total_mass
233 $ write(144,*)'Time ',time,'Mass ',total_mass
234 return
235 end
236c-----------------------------------------------------------------------
237! That was a bunch of diagostic routines
238c-----------------------------------------------------------------------
239 subroutine torque_calc_cmt(scale,x0,ifdout,iftout)
240c
241c Compute torque about point x0
242c
243c Scale is a user-supplied multiplier so that results may be
244c scaled to any convenient non-dimensionalization.
245c
246c
247 INCLUDE 'SIZE'
248 INCLUDE 'TOTAL'
249
250 common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
251 $ , scale_vf(3)
252
253c
254 real x0(3),w1(0:maxobj)
255 logical ifdout,iftout
256c
257 common /scrns/ sij (lx1*ly1*lz1*6*lelt)
258 common /scrcg/ pm1 (lx1,ly1,lz1,lelt)
259 common /scrsf/ xm0(lx1,ly1,lz1,lelt)
260 $, ym0(lx1,ly1,lz1,lelt)
261 $, zm0(lx1,ly1,lz1,lelt)
262c
263 parameter (lr=lx1*ly1*lz1)
264 common /scruz/ ur(lr),us(lr),ut(lr)
265 $ , vr(lr),vs(lr),vt(lr)
266 $ , wr(lr),ws(lr),wt(lr)
267c
268 common /ICPVARS/ pvars(lx1,ly1,lz1,7)
269 real pvars
270c
271 n = lx1*ly1*lz1*nelv
272 ntot1 = lx1*ly1*lz1
273c
274c call mappr(pm1,pr,xm0,ym0) ! map pressure onto Mesh 1
275c IN CMT-Nek Pressure is defined on the velocity grid
276 do ie=1,nelv
277 call copy(pm1(1,1,1,ie),pr(1,1,1,ie),ntot1)
278 enddo
279c
280c Add mean_pressure_gradient.X to p:
281
282c In CMT-Nek we do not fix the volume flow rate
283c if (param(55).ne.0) then
284c dpdx_mean = -scale_vf(1)
285c dpdy_mean = -scale_vf(2)
286c dpdz_mean = -scale_vf(3)
287c endif
288
289c call add2s2(pm1,xm1,dpdx_mean,n) ! Doesn't work if object is cut by
290c call add2s2(pm1,ym1,dpdy_mean,n) ! periodicboundary. In this case,
291c call add2s2(pm1,zm1,dpdz_mean,n) ! set ._mean=0 and compensate in
292c
293c Compute sij
294c
295 nij = 3
296 if (if3d.or.ifaxis) nij=6
297c we store primivtive variables only on dealiased grid. We will modify
298c comp_sij such that primvars are stored on the grid element
299c by element in a scratch array.
300 call comp_sij_cmt(sij,nij,ur,us,ut,vr,vs,vt,wr,ws,wt)
301c
302c
303c Fill up viscous array w/ default
304 param(2) = 0.0
305c
306 if (istep.lt.1) call cfill(vdiff,param(2),n)
307c
308 call cadd2(xm0,xm1,-x0(1),n)
309 call cadd2(ym0,ym1,-x0(2),n)
310 call cadd2(zm0,zm1,-x0(3),n)
311c
312 x1min=glmin(xm0(1,1,1,1),n)
313 x2min=glmin(ym0(1,1,1,1),n)
314 x3min=glmin(zm0(1,1,1,1),n)
315c
316 x1max=glmax(xm0(1,1,1,1),n)
317 x2max=glmax(ym0(1,1,1,1),n)
318 x3max=glmax(zm0(1,1,1,1),n)
319c
320 do i=0,maxobj
321 dragpx(i) = 0 ! BIG CODE :}
322 dragvx(i) = 0
323 dragx (i) = 0
324 dragpy(i) = 0
325 dragvy(i) = 0
326 dragy (i) = 0
327 dragpz(i) = 0
328 dragvz(i) = 0
329 dragz (i) = 0
330 torqpx(i) = 0
331 torqvx(i) = 0
332 torqx (i) = 0
333 torqpy(i) = 0
334 torqvy(i) = 0
335 torqy (i) = 0
336 torqpz(i) = 0
337 torqvz(i) = 0
338 torqz (i) = 0
339 enddo
340c
341c
342 nobj = 0
343 do ii=1,nhis
344 if (hcode(10,ii).EQ.'I') then
345 iobj = lochis(1,ii)
346 memtot = nmember(iobj)
347 nobj = max(iobj,nobj)
348c
349 if (hcode(1,ii).ne.' ' .or. hcode(2,ii).ne.' ' .or.
350 $ hcode(3,ii).ne.' ' ) then
351 ifield = 1
352c
353c Compute drag for this object
354c
355 do mem=1,memtot
356 ieg = object(iobj,mem,1)
357 ifc = object(iobj,mem,2)
358 if (gllnid(ieg).eq.nid) then ! this processor has a contribution
359 ie = gllel(ieg)
360 call drgtrq(dgtq,xm0,ym0,zm0,sij,pm1,vdiff,ifc,ie)
361c
362 call cmult(dgtq,scale,12)
363c
364 dragpx(iobj) = dragpx(iobj) + dgtq(1,1) ! pressure
365 dragpy(iobj) = dragpy(iobj) + dgtq(2,1)
366 dragpz(iobj) = dragpz(iobj) + dgtq(3,1)
367c
368 dragvx(iobj) = dragvx(iobj) + dgtq(1,2) ! viscous
369 dragvy(iobj) = dragvy(iobj) + dgtq(2,2)
370 dragvz(iobj) = dragvz(iobj) + dgtq(3,2)
371c
372 torqpx(iobj) = torqpx(iobj) + dgtq(1,3) ! pressure
373 torqpy(iobj) = torqpy(iobj) + dgtq(2,3)
374 torqpz(iobj) = torqpz(iobj) + dgtq(3,3)
375c
376 torqvx(iobj) = torqvx(iobj) + dgtq(1,4) ! viscous
377 torqvy(iobj) = torqvy(iobj) + dgtq(2,4)
378 torqvz(iobj) = torqvz(iobj) + dgtq(3,4)
379c
380 endif
381 enddo
382 endif
383 endif
384 enddo
385c
386c Sum contributions from all processors
387c
388 call gop(dragpx,w1,'+ ',maxobj+1)
389 call gop(dragpy,w1,'+ ',maxobj+1)
390 call gop(dragpz,w1,'+ ',maxobj+1)
391 call gop(dragvx,w1,'+ ',maxobj+1)
392 call gop(dragvy,w1,'+ ',maxobj+1)
393 call gop(dragvz,w1,'+ ',maxobj+1)
394c
395 call gop(torqpx,w1,'+ ',maxobj+1)
396 call gop(torqpy,w1,'+ ',maxobj+1)
397 call gop(torqpz,w1,'+ ',maxobj+1)
398 call gop(torqvx,w1,'+ ',maxobj+1)
399 call gop(torqvy,w1,'+ ',maxobj+1)
400 call gop(torqvz,w1,'+ ',maxobj+1)
401c
402 nobj = iglmax(nobj,1)
403c
404 do i=1,nobj
405 dragx(i) = dragpx(i) + dragvx(i)
406 dragy(i) = dragpy(i) + dragvy(i)
407 dragz(i) = dragpz(i) + dragvz(i)
408c
409 torqx(i) = torqpx(i) + torqvx(i)
410 torqy(i) = torqpy(i) + torqvy(i)
411 torqz(i) = torqpz(i) + torqvz(i)
412c
413 dragpx(0) = dragpx (0) + dragpx (i)
414 dragvx(0) = dragvx (0) + dragvx (i)
415 dragx (0) = dragx (0) + dragx (i)
416c
417 dragpy(0) = dragpy (0) + dragpy (i)
418 dragvy(0) = dragvy (0) + dragvy (i)
419 dragy (0) = dragy (0) + dragy (i)
420c
421 dragpz(0) = dragpz (0) + dragpz (i)
422 dragvz(0) = dragvz (0) + dragvz (i)
423 dragz (0) = dragz (0) + dragz (i)
424c
425 torqpx(0) = torqpx (0) + torqpx (i)
426 torqvx(0) = torqvx (0) + torqvx (i)
427 torqx (0) = torqx (0) + torqx (i)
428c
429 torqpy(0) = torqpy (0) + torqpy (i)
430 torqvy(0) = torqvy (0) + torqvy (i)
431 torqy (0) = torqy (0) + torqy (i)
432c
433 torqpz(0) = torqpz (0) + torqpz (i)
434 torqvz(0) = torqvz (0) + torqvz (i)
435 torqz (0) = torqz (0) + torqz (i)
436c
437 enddo
438c
439 i0 = 0
440 if (nobj.le.1) i0 = 1 ! one output for single-object case
441c
442 do i=i0,nobj
443 if (nio.eq.0) then
444 if (if3d.or.ifaxis) then
445 if (ifdout) then
446 write(6,6) istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
447 write(6,6) istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
448 write(6,6) istep,time,dragz(i),dragpz(i),dragvz(i),i,'dragz'
449 endif
450 if (iftout) then
451 write(6,6) istep,time,torqx(i),torqpx(i),torqvx(i),i,'torqx'
452 write(6,6) istep,time,torqy(i),torqpy(i),torqvy(i),i,'torqy'
453 write(6,6) istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
454 endif
455 else
456 if (ifdout) then
457 write(6,6) istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
458 write(6,6) istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
459 endif
460 if (iftout) then
461 write(6,6) istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
462 endif
463 endif
464 endif
465 6 format(i8,1p4e19.11,2x,i5,a5)
466 enddo
467c
468 return
469 end
470c-----------------------------------------------------------------------
471 subroutine comp_sij_cmt(sij,nij,ur,us,ut,vr,vs,vt,wr,ws,wt)
472c du_i du_j
473c Compute the stress tensor S_ij := ---- + ----
474c du_j du_i
475c
476 include 'SIZE'
477 include 'TOTAL'
478c
479 integer e
480c
481 real sij(lx1*ly1*lz1,nij,lelt)
482 real ur (1) , us (1) , ut (1)
483 $ , vr (1) , vs (1) , vt (1)
484 $ , wr (1) , ws (1) , wt (1)
485
486 real j ! Inverse Jacobian
487
488 common /ICPVARS/ pvars(lx1,ly1,lz1,7)
489 real pvars
490
491 n = lx1-1 ! Polynomial degree
492 nxyz = lx1*ly1*lz1
493
494 if (if3d) then ! 3D CASE
495 do e=1,nelv
496 call copy(pvars(1,1,1,1),vx(1,1,1,e),nxyz)
497 call copy(pvars(1,1,1,2),vy(1,1,1,e),nxyz)
498 call copy(pvars(1,1,1,3),vz(1,1,1,e),nxyz)
499 call local_grad3(ur,us,ut,pvars,N,1,dxm1,dxtm1)
500 call local_grad3(vr,vs,vt,pvars,N,2,dxm1,dxtm1)
501 call local_grad3(wr,ws,wt,pvars,N,3,dxm1,dxtm1)
502
503 do i=1,nxyz
504
505 j = jacmi(i,e)
506
507 sij(i,1,e) = j* ! du/dx + du/dx
508 $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e))
509
510 sij(i,2,e) = j* ! dv/dy + dv/dy
511 $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e))
512
513 sij(i,3,e) = j* ! dw/dz + dw/dz
514 $ 2*(wr(i)*rzm1(i,1,1,e)+ws(i)*szm1(i,1,1,e)+wt(i)*tzm1(i,1,1,e))
515
516 sij(i,4,e) = j* ! du/dy + dv/dx
517 $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e) +
518 $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e) )
519
520 sij(i,5,e) = j* ! dv/dz + dw/dy
521 $ (wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e) +
522 $ vr(i)*rzm1(i,1,1,e)+vs(i)*szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e) )
523
524 sij(i,6,e) = j* ! du/dz + dw/dx
525 $ (ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e) +
526 $ wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e) )
527
528 enddo
529 enddo
530
531 elseif (ifaxis) then ! AXISYMMETRIC CASE
532
533c
534c Notation: ( 2 x Acheson, p. 353)
535c Cylindrical
536c Nek5k Coordinates
537c
538c x z
539c y r
540c z theta
541c
542
543 do e=1,nelv
544 call copy(pvars(1,1,1,1),vx(1,1,1,e),nxyz)
545 call copy(pvars(1,1,1,2),vy(1,1,1,e),nxyz)
546 call copy(pvars(1,1,1,3),vz(1,1,1,e),nxyz)
547 call setaxdy ( ifrzer(e) ) ! change dytm1 if on-axis
548 call local_grad2(ur,us,pvars,N,1,dxm1,dytm1)
549 call local_grad2(vr,vs,pvars,N,2,dxm1,dytm1)
550 call local_grad2(wr,ws,pvars,N,3,dxm1,dytm1)
551
552 do i=1,nxyz
553 j = jacmi(i,e)
554 r = ym1(i,1,1,e) ! Cyl. Coord:
555
556 sij(i,1,e) = j* ! du/dx + du/dx ! e_zz
557 $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
558
559 sij(i,2,e) = j* ! dv/dy + dv/dy ! e_rr
560 $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
561
562 if (r.gt.0) then ! e_@@
563 sij(i,3,e) = pvars(i,1,1,1)/r ! v / r
564 else
565 sij(i,3,e) = j* ! L'Hopital's rule: e_@@ = dv/dr
566 $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
567 endif
568
569 sij(i,4,e) = j* ! du/dy + dv/dx ! e_zr
570 $ ( ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
571 $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
572
573 if (yyyr.gt.0) then ! e_r@
574 sij(i,5,e) = j* ! dw/dy
575 $ ( wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e) )
576 $ - pvars(i,1,1,3) / r
577 else
578 sij(i,5,e) = 0
579 endif
580
581 sij(i,6,e) = j* ! dw/dx ! e_@z
582 $ ( wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e) )
583
584 enddo
585 enddo
586
587 else ! 2D CASE
588
589 do e=1,nelv
590 call copy(pvars(1,1,1,1),vx(1,1,1,e),nxyz)
591 call copy(pvars(1,1,1,2),vy(1,1,1,e),nxyz)
592 call local_grad2(ur,us,pvars,N,1,dxm1,dxtm1)
593 call local_grad2(vr,vs,pvars,N,2,dxm1,dxtm1)
594
595 do i=1,nxyz
596 j = jacmi(i,e)
597
598 sij(i,1,e) = j* ! du/dx + du/dx
599 $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
600
601 sij(i,2,e) = j* ! dv/dy + dv/dy
602 $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
603
604 sij(i,3,e) = j* ! du/dy + dv/dx
605 $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
606 $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
607
608 enddo
609 enddo
610 endif
611 return
612 end
613c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.