source: CIVL/examples/fortran/nek5000/core/cmt/intpdiff.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: 8.8 KB
Line 
1C> @file intpdiff.f interpolation and differentiation routines not already provided
2C> by nek5000
3!--------------------------------------------------------------------
4! JH061914 propose name change to intpdiff since that is what is in
5! here.
6! JH081816 went to local_grad. redimensioned gradu. wish I could use
7! gradm11, but stride
8!--------------------------------------------------------------------
9
10 subroutine compute_gradients(e)
11 include 'SIZE'
12 include 'INPUT'
13 include 'DXYZ'
14 include 'GEOM'
15 include 'SOLN'
16 include 'CMTDATA'
17 parameter (ldd=lxd*lyd*lzd)
18 common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ud(ldd),tu(ldd)
19
20 integer eq,e
21
22! ! Compute d/dx, d/dy and d/dz of all the cons vars
23
24 nxy1 = lx1*ly1
25 nyz1 = ly1*lz1
26 nxyz1 = lx1*ly1*lz1
27 m0 = lx1-1
28
29 do eq=1,toteq
30 call invcol3(ud,u(1,1,1,eq,e),phig(1,1,1,e),nxyz1)
31
32 if (if3d) then
33 call local_grad3(ur,us,ut,ud,m0,1,dxm1,dxtm1)
34 do i=1,nxyz1
35 gradu(i,eq,1) = jacmi(i,e)*(rxm1(i,1,1,e)*ur(i)+
36 > sxm1(i,1,1,e)*us(i)+
37 > txm1(i,1,1,e)*ut(i))
38 enddo
39 do i=1,nxyz1
40 gradu(i,eq,2) = jacmi(i,e)*(rym1(i,1,1,e)*ur(i)+
41 > sym1(i,1,1,e)*us(i)+
42 > tym1(i,1,1,e)*ut(i))
43 enddo
44 do i=1,nxyz1
45 gradu(i,eq,3) = jacmi(i,e)*(rzm1(i,1,1,e)*ur(i)+
46 > szm1(i,1,1,e)*us(i)+
47 > tzm1(i,1,1,e)*ut(i))
48 enddo
49
50 else
51
52 call local_grad2(ur,us ,ud,m0,1,dxm1,dxtm1)
53 do i=1,nxyz1
54 gradu(i,eq,1) = jacmi(i,e)*(rxm1(i,1,1,e)*ur(i)+
55 > sxm1(i,1,1,e)*us(i))
56 enddo
57 do i=1,nxyz1
58 gradu(i,eq,2) = jacmi(i,e)*(rym1(i,1,1,e)*ur(i)+
59 > sym1(i,1,1,e)*us(i))
60 enddo
61
62 endif
63
64 enddo ! equation loop
65
66 return
67 end
68
69!-----------------------------------------------------------------------
70
71 subroutine set_dealias_face
72
73!-----------------------------------------------------------------------
74! JH111815 needed for face Jacobian and surface integrals
75!-----------------------------------------------------------------------
76
77 include 'SIZE'
78 include 'INPUT' ! for if3d
79 include 'GEOM' ! for ifgeom
80 include 'TSTEP' ! for istep
81 include 'WZ' ! for wxm1
82 include 'DG' ! for facewz
83
84 integer ilstep
85 save ilstep
86 data ilstep /-1/
87
88 if (.not.ifgeom.and.ilstep.gt.1) return ! already computed
89 if (ifgeom.and.ilstep.eq.istep) return ! already computed
90 ilstep = istep
91
92 call zwgl(zptf,wgtf,lxd)
93
94 if (if3d) then
95 k=0
96 do j=1,ly1
97 do i=1,lx1
98 k=k+1
99 wghtc(k)=wxm1(i)*wzm1(j)
100 enddo
101 enddo
102 k=0
103 do j=1,lyd
104 do i=1,lxd
105 k=k+1
106 wghtf(k)=wgtf(i)*wgtf(j)
107 enddo
108 enddo
109 else
110 call copy(wghtc,wxm1,lx1)
111 call copy(wghtf,wgtf,lxd)
112 endif
113
114 return
115 end
116!-----------------------------------------------------------------------
117
118 subroutine set_alias_rx(istp)
119! note that set_alias_rx will be called only when lxd = lx1
120 include 'SIZE'
121 include 'INPUT'
122 include 'GEOM'
123c include 'TSTEP' ! for istep
124
125 common /dealias1/ zd(lx1),wd(lx1)
126 integer e
127
128 integer ilsp
129 save ilsp
130 data ilsp /-1/
131
132 if (.not.ifgeom.and.ilsp.gt.1) return ! already computed
133 if (ifgeom.and.ilsp.eq.istp) return ! already computed
134 ilsp = istp
135 nxyz = lx1*ly1*lz1
136 call zwgll (zd,wd,lx1)
137
138 if (if3d)then
139 do e=1,nelv
140
141 call copy(rx(1,1,e),rxm1(1,1,1,e),nxyz)
142 call copy(rx(1,2,e),rym1(1,1,1,e),nxyz)
143 call copy(rx(1,3,e),rzm1(1,1,1,e),nxyz)
144 call copy(rx(1,4,e),sxm1(1,1,1,e),nxyz)
145 call copy(rx(1,5,e),sym1(1,1,1,e),nxyz)
146 call copy(rx(1,6,e),szm1(1,1,1,e),nxyz)
147 call copy(rx(1,7,e),txm1(1,1,1,e),nxyz)
148 call copy(rx(1,8,e),tym1(1,1,1,e),nxyz)
149 call copy(rx(1,9,e),tzm1(1,1,1,e),nxyz)
150
151 l = 0
152 do k=1,lz1
153 do j=1,ly1
154 do i=1,lx1
155 l = l+1
156 w = wd(i)*wd(j)*wd(k)
157 do ii=1,9
158 rx(l,ii,e) = w*rx(l,ii,e)
159 enddo
160 enddo
161 enddo
162 enddo
163 enddo
164 else
165
166 do e=1,nelv
167
168c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
169
170 call copy(rx(1,1,e),rxm1(1,1,1,e),nxyz)
171 call copy(rx(1,2,e),rym1(1,1,1,e),nxyz)
172 call copy(rx(1,3,e),sxm1(1,1,1,e),nxyz)
173 call copy(rx(1,4,e),sym1(1,1,1,e),nxyz)
174
175 l = 0
176 do j=1,ly1
177 do i=1,lx1
178 l = l+1
179 w = wd(i)*wd(j)
180 do ii=1,4
181 rx(l,ii,e) = w*rx(l,ii,e)
182 enddo
183 enddo
184 enddo
185 enddo
186
187 endif
188
189 return
190 end
191
192!-----------------------------------------------------------------------
193
194 subroutine gradm11_t(grad,uxyz,csgn,e) ! grad is incremented, not overwritten
195c source: . gradm1, from navier5.f
196c . gradm1, from navier5.f
197c Compute divergence^T of ux,uy,uz -- mesh 1 to mesh 1 (vel. to vel.)
198! single element, but references jacmi and the metrics
199
200 include 'SIZE'
201 include 'DXYZ'
202 include 'GEOM'
203 include 'INPUT'
204
205 parameter (lxyz=lx1*ly1*lz1)
206 real grad(lxyz),uxyz(lxyz,ldim)
207
208 common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz),ud(lxyz),tmp(lxyz)
209 real ur,us,ut,tmp
210
211 integer e
212
213 nxyz = lx1*ly1*lz1
214 call rzero(ud,nxyz)
215
216 N = lx1-1
217 if (if3d) then
218
219 do i=1,lxyz
220 ur(i) = jacmi(i,e)*(uxyz(i,1)*rxm1(i,1,1,e)
221 > + uxyz(i,2)*rym1(i,1,1,e)
222 > + uxyz(i,3)*rzm1(i,1,1,e) )
223 us(i) = jacmi(i,e)*(uxyz(i,1)*sxm1(i,1,1,e)
224 > + uxyz(i,2)*sym1(i,1,1,e)
225 > + uxyz(i,3)*szm1(i,1,1,e) )
226 ut(i) = jacmi(i,e)*(uxyz(i,1)*txm1(i,1,1,e)
227 > + uxyz(i,2)*tym1(i,1,1,e)
228 > + uxyz(i,3)*tzm1(i,1,1,e) )
229 enddo
230 call local_grad3_t(ud,ur,us,ut,N,1,dxm1,dxtm1,tmp)
231 else
232 do i=1,lxyz
233 ur(i) =jacmi(i,e)*(uxyz(i,1)*rxm1(i,1,1,e)
234 > + uxyz(i,2)*rym1(i,1,1,e) )
235 us(i) =jacmi(i,e)*(uxyz(i,1)*sxm1(i,1,1,e)
236 > + uxyz(i,2)*sym1(i,1,1,e) )
237 enddo
238 call local_grad2_t(ud,ur,us,N,1,dxm1,dxtm1,tmp)
239 endif
240 call cmult(ud,csgn,nxyz)
241 call add2(grad,ud,nxyz)
242
243 return
244 end
245
246!-----------------------------------------------------------------------
247
248 subroutine gradm1_t(u,ux,uy,uz)
249! JH082516 torn bleeding from Lu's dgf3.f. someday you are going to have
250! to vectorize cmt-nek properly.
251c source: . gradm1, from navier5.f
252c . gradm1, from navier5.f
253c
254c Compute divergence of ux,uy,uz -- mesh 1 to mesh 1 (vel. to vel.)
255c
256 include 'SIZE'
257 include 'DXYZ'
258 include 'GEOM'
259 include 'INPUT'
260c
261 parameter (lxyz=lx1*ly1*lz1)
262 real ux(lxyz,*),uy(lxyz,*),uz(lxyz,*),u(lxyz,*)
263
264 common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz),tmp(lxyz)
265 real ur,us,ut,tmp
266
267 integer e
268
269 nxyz = lx1*ly1*lz1
270 ntot = nxyz*nelt
271
272 N = lx1-1
273 do e=1,nelt
274 if (if3d) then
275
276 do i=1,lxyz
277 ur(i) = jacmi(i,e)*(ux(i,e)*rxm1(i,1,1,e)
278 $ + uy(i,e)*rym1(i,1,1,e)
279 $ + uz(i,e)*rzm1(i,1,1,e) )
280 us(i) = jacmi(i,e)*(ux(i,e)*sxm1(i,1,1,e)
281 $ + uy(i,e)*sym1(i,1,1,e)
282 $ + uz(i,e)*szm1(i,1,1,e) )
283 ut(i) = jacmi(i,e)*(ux(i,e)*txm1(i,1,1,e)
284 $ + uy(i,e)*tym1(i,1,1,e)
285 $ + uz(i,e)*tzm1(i,1,1,e) )
286 enddo
287 call local_grad3_t(u,ur,us,ut,N,e,dxm1,dxtm1,tmp)
288 else
289 do i=1,lxyz
290 ur(i) =jacmi(i,e)*(ux(i,e)*rxm1(i,1,1,e)
291 $ + uy(i,e)*rym1(i,1,1,e) )
292 us(i) =jacmi(i,e)*(ux(i,e)*sxm1(i,1,1,e)
293 $ + uy(i,e)*sym1(i,1,1,e) )
294 enddo
295 call local_grad2_t(u,ur,us,N,e,dxm1,dxtm1,tmp)
296 endif
297 enddo
298c
299 return
300 end
Note: See TracBrowser for help on using the repository browser.