source: CIVL/examples/fortran/nek5000/core/cmt/step.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: 11.1 KB
Line 
1C> @file step.f time stepping and mesh spacing routines
2 subroutine setdtcmt
3!--------------------------------------------------------------
4! JH072914 now summed instead of maximized for compressible flow
5! JH082015 why was I doing that again?
6! someday compute new timestep based on CFL-condition. follow
7! setdtc in nek/subs1.f very carefully.
8! JH091616 now with diffusion number for naively explicit visc
9! JH091616 consider migrating to compute_cfl
10! JH120116 Why aren't we including total again?
11!--------------------------------------------------------------
12 include 'SIZE'
13 include 'GEOM'
14 include 'MVGEOM'
15 include 'MASS'
16 include 'TSTEP'
17 include 'INPUT'
18 include 'SOLN'
19 include 'CMTDATA'
20
21!--------------------------------------------------------------
22! YOU REALLY PROBABLY WANT YOUR OWN SCRATCH SPACE FOR THIS
23 common /scrsf/ utmp(lx1,ly1,lz1,lelt) ! mind if I borrow this?
24 $ , vtmp(lx1,ly1,lz1,lelt) ! as long as the mesh
25 $ , wtmp(lx1,ly1,lz1,lelt) ! doesn't move
26! YOU REALLY PROBABLY WANT YOUR OWN SCRATCH SPACE FOR THAT
27!--------------------------------------------------------------
28 common /udxmax/ umax
29 real strof
30 data strof /1.0e-8/
31
32 dt=abs(param(12))
33
34 NTOT = lx1*ly1*lz1*NELV
35 do i=1,ntot
36 utmp(i,1,1,1) = abs(vx(i,1,1,1))+csound(i,1,1,1)
37 vtmp(i,1,1,1) = abs(vy(i,1,1,1))+csound(i,1,1,1)
38 wtmp(i,1,1,1) = abs(vz(i,1,1,1))+csound(i,1,1,1)
39 enddo
40
41! JH091118 DefaultParameters means we don't have direct control over
42! this variable anymore. If it's lower than its default value,
43! we trust it to set the time step
44! if (ctarg .gt.0.0) then
45 if (ctarg .lt.0.5) then
46 call compute_cfl (umax,utmp,vtmp,wtmp,1.0)
47 dt_cfl=ctarg/umax
48 call glsqinvcolmin(dt1,vdiff(1,1,1,1,imu ),gridh,ntot,ctarg)
49 call glsqinvcolmin(dt2,vdiff(1,1,1,1,iknd),gridh,ntot,ctarg)
50 call glsqinvcolmin(dt3,vdiff(1,1,1,1,inus),gridh,ntot,ctarg)
51 dt=min(dt_cfl,dt1,dt2,dt3)
52 if (dt .gt. 10.0) then
53 if (nio.eq.0) write(6,*) 'dt huge. crashing ',istep,stage,
54 > dt
55 call exitt
56 endif
57 else
58 dt = abs(param(12))
59 endif
60 dt_ptcle = dt
61#ifdef LPM
62 call lpm_set_dt(dt_ptcle) ! particle time step
63 dt=min(dt,dt_ptcle)
64#endif
65
66 if (timeio .gt. 0.0) then ! adjust dt for timeio.
67 zetime1=time_cmt
68 zetime2=time_cmt+dt
69 it1=zetime1/timeio
70 it2=zetime2/timeio
71 ita=it1
72 itb=ita+1
73 if (abs(zetime1-itb*timeio).le.strof) it1=itb
74 ita=it2
75 itb=ita+1
76 if (abs(zetime2-itb*timeio).le.strof) it2=itb
77 if (it2.gt.it1) then
78 ifoutfld=.true.
79 dt=(it2*timeio)-time_cmt
80 endif
81 endif
82
83 param(12)=-dt
84
85 call compute_cfl (courno,utmp,vtmp,wtmp,dt) ! sanity?
86
87! diffusion number based on viscosity.
88
89! call mindr(mdr,diffno2)
90 call glinvcol2max(diffno1,vdiff(1,1,1,1,imu), gridh,ntot,dt)
91 call glinvcol2max(diffno2,vdiff(1,1,1,1,iknd),gridh,ntot,dt)
92 call glinvcol2max(diffno3,vdiff(1,1,1,1,inus),gridh,ntot,dt)
93! diffno=max(diffno1,diffno2,diffno3)
94 time_cmt= time_cmt+dt
95 time = time_cmt
96 if (nio.eq.0) WRITE(6,100)ISTEP,TIME_CMT,DT,COURNO,
97 > diffno1,diffno2,diffno3
98 100 FORMAT('CMT ',I7,', t=',1pE14.7,', DT=',1pE14.7
99 $,', C=',1pE12.5,', Dmu,knd,art=',3(1pE11.4))
100
101 return
102 end
103
104 subroutine mindr(mdr,diffno)
105c
106c Find minimum distance between grid points
107c and multiply it by viscosity to get diffusion number
108c Probably need to do this for energy equation too...
109! JH091616 migrate to getdr 3d ASAP. Again, follow compute_cfl
110 include 'SIZE'
111 include 'TOTAL'
112 include 'CMTDATA'
113
114 real mdr,dr1
115 real x0,x1,x2,x3,x4,y0,y1,y2,y3,y4
116 integer e
117c
118 mdr=1.e5
119 if(if3d) then
120 write(6,*)'TODO:: mindr for 3D. getdr'
121 call exitt
122 else
123 diffno=0.0
124 do e=1,nelt
125 do iy=2,ly1-1
126 do ix=2,lx1-1
127 dtmp=1.0e5
128 x0 = xm1(ix ,iy ,1,e)
129 x1 = xm1(ix ,iy-1,1,e)
130 x2 = xm1(ix+1,iy ,1,e)
131 x3 = xm1(ix ,iy+1,1,e)
132 x4 = xm1(ix-1,iy ,1,e)
133 y0 = ym1(ix ,iy ,1,e)
134 y1 = ym1(ix ,iy-1,1,e)
135 y2 = ym1(ix+1,iy ,1,e)
136 y3 = ym1(ix ,iy+1,1,e)
137 y4 = ym1(ix-1,iy ,1,e)
138 dr1 = dist2(x0,y0,x1,y1)
139 if(dr1.lt.mdr) mdr=dr1
140 if(dr1.lt.dtmp) dtmp=dr1
141 dr1 = dist2(x0,y0,x2,y2)
142 if(dr1.lt.mdr) mdr=dr1
143 if(dr1.lt.dtmp) dtmp=dr1
144 dr1 = dist2(x0,y0,x3,y3)
145 if(dr1.lt.mdr) mdr=dr1
146 if(dr1.lt.dtmp) dtmp=dr1
147 dr1 = dist2(x0,y0,x4,y4)
148 if(dr1.lt.mdr) mdr=dr1
149 if(dr1.lt.dtmp) dtmp=dr1
150 diffno=max(diffno,dt*vdiff(ix,iy,1,e,imu)/dtmp/dtmp)
151 enddo
152 enddo
153 enddo
154 endif
155 diffno=glmax(diffno,1)
156 mdr = glmin(mdr,1)
157 if(mdr.ge.1.e5) write(6,*) 'wrong mdr'
158
159 return
160 end
161
162 real function dist2(x1,y1,x2,y2)
163 real x1,y1,x2,y2,dx,dy
164c
165 dx = x1-x2
166 dy = y1-y2
167 dist2 = sqrt(dx*dx+dy*dy)
168c
169 return
170 end
171
172!-----------------------------------------------------------------------
173
174 subroutine compute_grid_h(h,x,y,z)
175! Richard Pasquetti SEM "grid spacing h." good parallelogram/piped stuff
176 include 'SIZE'
177 include 'INPUT'
178 include 'GEOM'
179 real a(3), b(3), c(3), d(3)
180 real h(lx1,ly1,lz1,nelt) ! intent(out)
181 real x(lx1,ly1,lz1,nelt) ! intent(in)
182 real y(lx1,ly1,lz1,nelt) ! intent(in)
183 real z(lx1,ly1,lz1,nelt) ! intent(in)
184 integer e
185 integer icalld
186 data icalld /0/
187 save icalld
188
189 if (icalld .eq. 1) then
190 return
191 else
192 icalld=1
193 endif
194
195 do e=1,nelt
196 do iz=1,lz1
197 if (if3d) then
198 km1=iz-1
199 kp1=iz+1
200 izm=km1
201 if (km1 .lt. 1) izm=iz
202 izp=kp1
203 if (kp1 .gt. lz1) izp=iz
204 else
205 izm=iz
206 izp=iz
207 endif
208 do iy=1,ly1
209 jm1=iy-1
210 jp1=iy+1
211 iym=jm1
212 if (jm1 .lt. 1) iym=iy
213 iyp=jp1
214 if (jp1 .gt. ly1) iyp=iy
215 do ix=1,lx1
216 im1=ix-1
217 ip1=ix+1
218 ixm=im1
219 if (im1 .lt. 1) ixm=ix
220 ixp=ip1
221 if (ip1 .gt. lx1) ixp=ix
222 x1 = x(ixm,iy ,iz ,e)
223 x2 = x(ixp,iy ,iz ,e)
224 x3 = x(ix ,iym,iz ,e)
225 x4 = x(ix ,iyp,iz ,e)
226 x5 = x(ix ,iy ,izm,e)
227 x6 = x(ix ,iy ,izp,e)
228 y1 = y(ixm,iy ,iz ,e)
229 y2 = y(ixp,iy ,iz ,e)
230 y3 = y(ix ,iym,iz ,e)
231 y4 = y(ix ,iyp,iz ,e)
232 y5 = y(ix ,iy ,izm,e)
233 y6 = y(ix ,iy ,izp,e)
234 z1 = z(ixm,iy ,iz ,e)
235 z2 = z(ixp,iy ,iz ,e)
236 z3 = z(ix ,iym,iz ,e)
237 z4 = z(ix ,iyp,iz ,e)
238 z5 = z(ix ,iy ,izm,e)
239 z6 = z(ix ,iy ,izp,e)
240 a(1)=x2-x1
241 a(2)=y2-y1
242 a(3)=z2-z1
243 b(1)=x4-x3
244 b(2)=y4-y3
245 b(3)=z4-z3
246 c(1)=x6-x5
247 c(2)=y6-y5
248 c(3)=z6-z5
249 if (if3d) then
250 fact=0.125 ! h doesn't reach into corners of neighboring elements
251 if (ixp.eq.ix.or.ixm.eq.ix) fact=2.0*fact
252 if (iym.eq.iy.or.iyp.eq.iy) fact=2.0*fact
253 if (izm.eq.iz.or.izp.eq.iz) fact=2.0*fact
254 call cross(d,a,b)
255 h(ix,iy,iz,e)=fact*dot(c,d,3)
256 h(ix,iy,iz,e)=abs(h(ix,iy,iz,e))**(1.0/3.0)
257 else
258 fact=0.25
259 if (ixp.eq.ix.or.ixm.eq.ix) fact=2.0*fact
260 if (iym.eq.iy.or.iyp.eq.iy) fact=2.0*fact
261 h(ix,iy,iz,e)=sqrt(fact*abs(a(1)*b(2)-a(2)*b(1)))
262 endif
263 enddo
264 enddo
265 enddo
266 enddo
267
268 return
269 end
270
271!-----------------------------------------------------------------------
272
273 subroutine glinvcol2max(col2m,a,b,n,s)
274 real col2m
275 real s
276 real a(*),b(*)
277 tmp=0.0
278 do i=1,n
279 tmp=max(tmp,abs(s*a(i)/b(i)/b(i)))
280 enddo
281 col2m=glamax(tmp,1)
282 return
283 end
284
285!-----------------------------------------------------------------------
286
287 subroutine glsqinvcolmin(col2m,a,b,n,s)
288 real col2m
289 real s
290 real a(*),b(*)
291 tmp=1.0e36
292 do i=1,n
293 if (a(i).gt.1.0e-36) tmp=min(tmp,abs(s*b(i)**2/a(i)))
294 enddo
295 col2m=glamin(tmp,1)
296 return
297 end
298
299!-----------------------------------------------------------------------
300
301 subroutine compute_mesh_h(h,x,y,z)
302! Zingan's DGFEM formula: h=minimum distance between vertices divided by
303! polynomial order
304 include 'SIZE'
305 include 'INPUT'
306 real h(nelt) ! intent(out)
307 real x(lx1,ly1,lz1,nelt) ! intent(in)
308 real y(lx1,ly1,lz1,nelt) ! intent(in)
309 real z(lx1,ly1,lz1,nelt) ! intent(in)
310 real xcrn(8),ycrn(8),zcrn(8)
311 integer e
312 integer icalld
313 data icalld /0/
314 save icalld
315
316 if (icalld .eq. 1) then
317 return
318 else
319 icalld=1
320 endif
321
322 ncrn=2**ldim
323 rp=1.0/((lx1-1))
324
325 do e=1,nelt
326 call rzero(zcrn,8)
327 k1=1
328 k2=lz1
329 j1=1
330 j2=ly1
331 i1=1
332 i2=lx1
333 xcrn(1) = x(i1,j1,k1,e)
334 xcrn(2) = x(i2,j1,k1,e)
335 xcrn(3) = x(i1,j2,k1,e)
336 xcrn(4) = x(i2,j2,k1,e)
337 ycrn(1) = y(i1,j1,k1,e)
338 ycrn(2) = y(i2,j1,k1,e)
339 ycrn(3) = y(i1,j2,k1,e)
340 ycrn(4) = y(i2,j2,k1,e)
341 if (if3d) then
342 xcrn(5) = x(i1,j1,k2,e)
343 xcrn(6) = x(i2,j1,k2,e)
344 xcrn(7) = x(i1,j2,k2,e)
345 xcrn(8) = x(i2,j2,k2,e)
346 ycrn(5) = y(i1,j1,k2,e)
347 ycrn(6) = y(i2,j1,k2,e)
348 ycrn(7) = y(i1,j2,k2,e)
349 ycrn(8) = y(i2,j2,k2,e)
350 zcrn(1) = z(i1,j1,k1,e)
351 zcrn(2) = z(i2,j1,k1,e)
352 zcrn(3) = z(i1,j2,k1,e)
353 zcrn(4) = z(i2,j2,k1,e)
354 zcrn(5) = z(i1,j1,k2,e)
355 zcrn(6) = z(i2,j1,k2,e)
356 zcrn(7) = z(i1,j2,k2,e)
357 zcrn(8) = z(i2,j2,k2,e)
358 endif
359 dist=1.0e36
360 do ic1=1,ncrn
361 do ic2=1,ncrn
362 if (ic2 .ne. ic1) then
363 dtmp=(xcrn(ic2)-xcrn(ic1))**2+
364 > (ycrn(ic2)-ycrn(ic1))**2+
365 > (zcrn(ic2)-zcrn(ic1))**2
366 dist=min(dist,sqrt(dtmp))
367 endif
368 enddo
369 enddo
370 h(e)=dist*rp
371 enddo
372
373 return
374 end
Note: See TracBrowser for help on using the repository browser.