CG  Version 25
insdt.h
Go to the documentation of this file.
1 c ===============================================================================
2 c This file is included by
3 c insdtINS.bf
4 c insdtKE.bf
5 c insdtSPAL.bf
6 c insdtVP.bf
7 c ==============================================================================
8 
9 c These next include files will define the macros that will define the difference approximations
10 c The actual macro is called below
11 #Include "defineDiffOrder2f.h"
12 #Include "defineDiffOrder4f.h"
13 
14 #Include "commonMacros.h"
15 
16 #beginMacro loopse1(e1)
17 if( useWhereMask.ne.0 )then
18  do i3=n3a,n3b
19  do i2=n2a,n2b
20  do i1=n1a,n1b
21  if( mask(i1,i2,i3).gt.0 )then
22  e1
23  end if
24  end do
25  end do
26  end do
27 else
28  do i3=n3a,n3b
29  do i2=n2a,n2b
30  do i1=n1a,n1b
31  e1
32  end do
33  end do
34  end do
35 end if
36 #endMacro
37 #beginMacro loopse2(e1,e2)
38 if( useWhereMask.ne.0 )then
39  do i3=n3a,n3b
40  do i2=n2a,n2b
41  do i1=n1a,n1b
42  if( mask(i1,i2,i3).gt.0 )then
43  e1
44  e2
45  end if
46  end do
47  end do
48  end do
49 else
50  do i3=n3a,n3b
51  do i2=n2a,n2b
52  do i1=n1a,n1b
53  e1
54  e2
55  end do
56  end do
57  end do
58 end if
59 #endMacro
60 
61 #beginMacro loopse3(e1,e2,e3)
62 if( useWhereMask.ne.0 )then
63  do i3=n3a,n3b
64  do i2=n2a,n2b
65  do i1=n1a,n1b
66  if( mask(i1,i2,i3).gt.0 )then
67  e1
68  e2
69  e3
70  end if
71  end do
72  end do
73  end do
74 else
75  do i3=n3a,n3b
76  do i2=n2a,n2b
77  do i1=n1a,n1b
78  e1
79  e2
80  e3
81  end do
82  end do
83  end do
84 end if
85 #endMacro
86 
87 #beginMacro loopse4(e1,e2,e3,e4)
88 if( useWhereMask.ne.0 )then
89  do i3=n3a,n3b
90  do i2=n2a,n2b
91  do i1=n1a,n1b
92  if( mask(i1,i2,i3).gt.0 )then
93  e1
94  e2
95  e3
96  e4
97  end if
98  end do
99  end do
100  end do
101 else
102  do i3=n3a,n3b
103  do i2=n2a,n2b
104  do i1=n1a,n1b
105  e1
106  e2
107  e3
108  e4
109  end do
110  end do
111  end do
112 end if
113 #endMacro
114 
115 #beginMacro loopse6(e1,e2,e3,e4,e5,e6)
116 if( useWhereMask.ne.0 )then
117  do i3=n3a,n3b
118  do i2=n2a,n2b
119  do i1=n1a,n1b
120  if( mask(i1,i2,i3).gt.0 )then
121  e1
122  e2
123  e3
124  e4
125  e5
126  e6
127  end if
128  end do
129  end do
130  end do
131 else
132  do i3=n3a,n3b
133  do i2=n2a,n2b
134  do i1=n1a,n1b
135  e1
136  e2
137  e3
138  e4
139  e5
140  e6
141  end do
142  end do
143  end do
144 end if
145 #endMacro
146 
147 
148 
149 
150 c Define the artificial diffusion coefficients
151 c gt should be R or C (gridType is Rectangular or Curvilinear)
152 c tb should be blank or SA (SA=Spalart-Allamras turbulence model)
153 #beginMacro defineArtificialDiffusionCoefficients(dim,gt,tb)
154  #If #dim == "2"
155  cdmz=admz ## gt ## tb(i1 ,i2 ,i3)
156  cdpz=admz ## gt ## tb(i1+1,i2 ,i3)
157  cdzm=adzm ## gt ## tb(i1 ,i2 ,i3)
158  cdzp=adzm ## gt ## tb(i1 ,i2+1,i3)
159  cdDiag=cdmz+cdpz+cdzm+cdzp
160  ! write(*,'(1x,''insdt:i1,i2,cdmz,cdzm='',2i3,2f9.3)') i1,i2,cdmz,cdzm
161  #Else
162  cdmzz=admzz ## gt ## tb(i1 ,i2 ,i3 )
163  cdpzz=admzz ## gt ## tb(i1+1,i2 ,i3 )
164  cdzmz=adzmz ## gt ## tb(i1 ,i2 ,i3 )
165  cdzpz=adzmz ## gt ## tb(i1 ,i2+1,i3 )
166  cdzzm=adzzm ## gt ## tb(i1 ,i2 ,i3 )
167  cdzzp=adzzm ## gt ## tb(i1 ,i2 ,i3+1)
168  cdDiag=cdmzz+cdpzz+cdzmz+cdzpz+cdzzm+cdzzp
169  #End
170 #endMacro
171 
172 c Define macros for the derivatives based on the dimension, order of accuracy and grid-type
173 #beginMacro defineDerivativeMacros(DIM,ORDER,GRIDTYPE)
174 
175 #defineMacro U(cc) u(i1,i2,i3,cc)
176 #defineMacro UU(cc) uu(i1,i2,i3,cc)
177 
178 #If #DIM == "2"
179  #If #ORDER == "2"
180  #If #GRIDTYPE == "rectangular"
181  #defineMacro UX(cc) ux22r(i1,i2,i3,cc)
182  #defineMacro UY(cc) uy22r(i1,i2,i3,cc)
183  #defineMacro UXX(cc) uxx22r(i1,i2,i3,cc)
184  #defineMacro UXY(cc) uxy22r(i1,i2,i3,cc)
185  #defineMacro UYY(cc) uyy22r(i1,i2,i3,cc)
186  #defineMacro ULAP(cc) ulaplacian22r(i1,i2,i3,cc)
187  #Else
188  #defineMacro UX(cc) ux22(i1,i2,i3,cc)
189  #defineMacro UY(cc) uy22(i1,i2,i3,cc)
190  #defineMacro UXX(cc) uxx22(i1,i2,i3,cc)
191  #defineMacro UXY(cc) uxy22(i1,i2,i3,cc)
192  #defineMacro UYY(cc) uyy22(i1,i2,i3,cc)
193  #defineMacro ULAP(cc) ulaplacian22(i1,i2,i3,cc)
194  #End
195  #Else
196  #If #GRIDTYPE == "rectangular"
197  #defineMacro UX(cc) ux42r(i1,i2,i3,cc)
198  #defineMacro UY(cc) uy42r(i1,i2,i3,cc)
199  #defineMacro UXX(cc) uxx42r(i1,i2,i3,cc)
200  #defineMacro UXY(cc) uxy42r(i1,i2,i3,cc)
201  #defineMacro UYY(cc) uyy42r(i1,i2,i3,cc)
202  #defineMacro ULAP(cc) ulaplacian42r(i1,i2,i3,cc)
203  #Else
204  #defineMacro UX(cc) ux42(i1,i2,i3,cc)
205  #defineMacro UY(cc) uy42(i1,i2,i3,cc)
206  #defineMacro UXX(cc) uxx42(i1,i2,i3,cc)
207  #defineMacro UXY(cc) uxy42(i1,i2,i3,cc)
208  #defineMacro UYY(cc) uyy42(i1,i2,i3,cc)
209  #defineMacro ULAP(cc) ulaplacian42(i1,i2,i3,cc)
210  #End
211  #End
212 #Else
213  #If #ORDER == "2"
214  #If #GRIDTYPE == "rectangular"
215  #defineMacro UX(cc) ux23r(i1,i2,i3,cc)
216  #defineMacro UY(cc) uy23r(i1,i2,i3,cc)
217  #defineMacro UZ(cc) uz23r(i1,i2,i3,cc)
218  #defineMacro UXX(cc) uxx23r(i1,i2,i3,cc)
219  #defineMacro UXY(cc) uxy23r(i1,i2,i3,cc)
220  #defineMacro UXZ(cc) uxz23r(i1,i2,i3,cc)
221  #defineMacro UYY(cc) uyy23r(i1,i2,i3,cc)
222  #defineMacro UYZ(cc) uyz23r(i1,i2,i3,cc)
223  #defineMacro UZZ(cc) uzz23r(i1,i2,i3,cc)
224  #defineMacro ULAP(cc) ulaplacian23r(i1,i2,i3,cc)
225  #Else
226  #defineMacro UX(cc) ux23(i1,i2,i3,cc)
227  #defineMacro UY(cc) uy23(i1,i2,i3,cc)
228  #defineMacro UZ(cc) uz23(i1,i2,i3,cc)
229  #defineMacro UXX(cc) uxx23(i1,i2,i3,cc)
230  #defineMacro UXY(cc) uxy23(i1,i2,i3,cc)
231  #defineMacro UXZ(cc) uxz23(i1,i2,i3,cc)
232  #defineMacro UYY(cc) uyy23(i1,i2,i3,cc)
233  #defineMacro UYZ(cc) uyz23(i1,i2,i3,cc)
234  #defineMacro UZZ(cc) uzz23(i1,i2,i3,cc)
235  #defineMacro ULAP(cc) ulaplacian23(i1,i2,i3,cc)
236  #End
237 
238  #Else
239 
240  #If #GRIDTYPE == "rectangular"
241  #defineMacro UX(cc) ux43r(i1,i2,i3,cc)
242  #defineMacro UY(cc) uy43r(i1,i2,i3,cc)
243  #defineMacro UZ(cc) uz43r(i1,i2,i3,cc)
244  #defineMacro UXX(cc) uxx43r(i1,i2,i3,cc)
245  #defineMacro UXY(cc) uxy43r(i1,i2,i3,cc)
246  #defineMacro UXZ(cc) uxz43r(i1,i2,i3,cc)
247  #defineMacro UYY(cc) uyy43r(i1,i2,i3,cc)
248  #defineMacro UYZ(cc) uyz43r(i1,i2,i3,cc)
249  #defineMacro UZZ(cc) uzz43r(i1,i2,i3,cc)
250  #defineMacro ULAP(cc) ulaplacian43r(i1,i2,i3,cc)
251  #Else
252  #defineMacro UX(cc) ux43(i1,i2,i3,cc)
253  #defineMacro UY(cc) uy43(i1,i2,i3,cc)
254  #defineMacro UZ(cc) uz43(i1,i2,i3,cc)
255  #defineMacro UXX(cc) uxx43(i1,i2,i3,cc)
256  #defineMacro UXY(cc) uxy43(i1,i2,i3,cc)
257  #defineMacro UXZ(cc) uxz43(i1,i2,i3,cc)
258  #defineMacro UYY(cc) uyy43(i1,i2,i3,cc)
259  #defineMacro UYZ(cc) uyz43(i1,i2,i3,cc)
260  #defineMacro UZZ(cc) uzz43(i1,i2,i3,cc)
261  #defineMacro ULAP(cc) ulaplacian43(i1,i2,i3,cc)
262  #End
263  #End
264 #End
265 #endMacro
266 
267 
268 c =============================================================
269 c Compute derivatives of u,v,w
270 c =============================================================
271 #beginMacro setupDerivatives(DIM)
272  u0x=UX(uc)
273  u0y=UY(uc)
274  v0x=UX(vc)
275  v0y=UY(vc)
276 #If #DIM == "3"
277  u0z=UZ(uc)
278  v0z=UZ(vc)
279  w0x=UX(wc)
280  w0y=UY(wc)
281  w0z=UZ(wc)
282 #End
283 #endMacro
284 
285 
286 
287 
288 
289 c***************************************************************
290 c Define the equations for EXPLICIT time stepping
291 c
292 c ORDER: 2,4
293 c DIM: 2,3
294 c GRIDTYPE: rectangular, curvilinear
295 c
296 c***************************************************************
297 #beginMacro fillEquations(DIM,ORDER,GRIDTYPE)
298 
299 defineDerivativeMacros(DIM,ORDER,GRIDTYPE)
300 
301 if( isAxisymmetric.eq.0 )then
302 
303  if( gridIsImplicit.eq.0 )then
304  ! explicit
305 
306  loopse1($buildEquations(EXPLICIT,NONE,DIM,ORDER,GRIDTYPE,NO))
307 
308  else ! gridIsImplicit
309  ! ***** implicit *******
310 
311  if( implicitOption .eq.computeImplicitTermsSeparately )then
312  loopse1($buildEquations(BOTH,NONE,DIM,ORDER,GRIDTYPE,NO))
313  else if( implicitOption.eq.doNotComputeImplicitTerms )then
314  loopse1($buildEquations(EXPLICIT_ONLY,NONE,DIM,ORDER,GRIDTYPE,NO))
315  else
316  write(*,*)'insdt: Unknown implicitOption=',implicitOption
317  stop 5
318  end if ! end implicitOption
319 
320  end if
321 
322 else if( isAxisymmetric.eq.1 )then
323 
324  #If (#DIM == "2")
325  ! **** axisymmetric case ****
326  if( gridIsImplicit.eq.0 )then
327  ! explicit
328 
329  loopse1($buildEquations(EXPLICIT,NONE,DIM,ORDER,GRIDTYPE,YES))
330 
331  else ! gridIsImplicit
332  ! ***** implicit *******
333 
334  if( implicitOption .eq.computeImplicitTermsSeparately )then
335  loopse1($buildEquations(BOTH,NONE,DIM,ORDER,GRIDTYPE,YES))
336  else if( implicitOption.eq.doNotComputeImplicitTerms )then
337  loopse1($buildEquations(EXPLICIT_ONLY,NONE,DIM,ORDER,GRIDTYPE,YES))
338  else
339  write(*,*)'insdt: Unknown implicitOption=',implicitOption
340  stop 5
341  end if ! end implicitOption
342 
343  end if
344  #End
345 
346 else
347  stop 88733
348 end if
349 
350 #endMacro
351 
352 
353 
354 c====================================================================================
355 c
356 c====================================================================================
357 #beginMacro assignEquations(DIM,ORDER)
358 if( gridType.eq.rectangular )then
359  fillEquations(DIM,ORDER,rectangular)
360 else if( gridType.eq.curvilinear )then
361  fillEquations(DIM,ORDER,curvilinear)
362 else
363  stop 77
364 end if
365 #endMacro
366 
367 
368 c======================================================================================
369 c Define the subroutine to compute du/dt
370 c
371 c======================================================================================
372 #beginMacro INSDT(NAME,DIM,ORDER)
373  subroutine NAME(nd,n1a,n1b,n2a,n2b,n3a,n3b,nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b,\
374  mask,xy,rsxy,radiusInverse, u,uu, ut,uti,gv,dw, bc, ipar, rpar, ierr )
375 c======================================================================
376 c Compute du/dt for the incompressible NS on rectangular grids
377 c OPTIMIZED version for rectangular grids.
378 c nd : number of space dimensions
379 c
380 c gv : gridVelocity for moving grids
381 c uu : for moving grids uu is a workspace to hold u-gv, otherwise uu==u
382 c dw : distance to the wall for some turbulence models
383 c======================================================================
384  implicit none
385  integer nd, n1a,n1b,n2a,n2b,n3a,n3b,nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b
386 
387  real u(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
388  real uu(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
389  real ut(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
390  real uti(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
391  real gv(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
392  real dw(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
393  real xy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
394  real rsxy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1,0:nd-1)
395  real radiusInverse(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
396 
397  integer mask(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
398  integer bc(0:1,0:2),ierr
399 
400  integer ipar(0:*)
401  real rpar(0:*)
402 
403  ! ---- local variables -----
404  integer c,i1,i2,i3,kd,kd3,orderOfAccuracy,gridIsMoving,useWhereMask
405  integer gridIsImplicit,implicitOption,implicitMethod,isAxisymmetric,use2ndOrderAD,use4thOrderAD
406  integer rc,pc,uc,vc,wc,sc,nc,kc,ec,tc,grid,m,advectPassiveScalar,vsc
407  real nu,dt,nuPassiveScalar,adcPassiveScalar
408  real dxi,dyi,dzi,dri,dsi,dti,dr2i,ds2i,dt2i
409  real ad21,ad22,ad41,ad42,cd22,cd42,adc
410  real ad21n,ad22n,ad41n,ad42n,cd22n,cd42n
411  real yy,ri
412 
413  integer gridType
414  integer rectangular,curvilinear
415  parameter( rectangular=0, curvilinear=1 )
416 
417  integer turbulenceModel,noTurbulenceModel
418  integer baldwinLomax,spalartAllmaras,kEpsilon,kOmega,largeEddySimulation
419  parameter (noTurbulenceModel=0,baldwinLomax=1,kEpsilon=2,kOmega=3,spalartAllmaras=4,largeEddySimulation=5 )
420 
421  integer pdeModel,standardModel,BoussinesqModel,viscoPlasticModel,twoPhaseFlowModel
422  parameter( standardModel=0,BoussinesqModel=1,viscoPlasticModel=2,twoPhaseFlowModel=3 )
423 
424  integer computeAllTerms,\
425  doNotComputeImplicitTerms,\
426  computeImplicitTermsSeparately,\
427  computeAllWithWeightedImplicit
428 
429  parameter( computeAllTerms=0,\
430  doNotComputeImplicitTerms=1,\
431  computeImplicitTermsSeparately=2,\
432  computeAllWithWeightedImplicit=3 )
433 
434  real rx,ry,rz,sx,sy,sz,tx,ty,tz
435  real dr(0:2), dx(0:2)
436 
437  ! for SPAL TM
438  real n0,n0x,n0y,n0z
439  real cb1, cb2, cv1, sigma, sigmai, kappa, cw1, cw2, cw3, cw3e6, cv1e3, cd0, cr0
440  real chi,chi3,fnu1,fnu2,s,r,g,fw,dKappaSq,nSqBydSq,dd
441  real nuT,nuTx,nuTy,nuTz,nuTd
442 
443  real u0,u0x,u0y,u0z
444  real v0,v0x,v0y,v0z
445  real w0,w0x,w0y,w0z
446  ! for k-epsilon
447  real k0,k0x,k0y,k0z, e0,e0x,e0y,e0z
448  real nuP,prod
449  real cMu,cEps1,cEps2,sigmaEpsI,sigmaKI
450 
451  ! for visco-plastic
452  ! real nuVP,etaVP,yieldStressVP,exponentVP,epsVP
453  ! real eDotNorm,exp0
454  ! real u0xx,u0xy,u0xz,u0yy,u0yz,u0zz
455  ! real v0xx,v0xy,v0xz,v0yy,v0yz,v0zz
456  ! real w0xx,w0xy,w0xz,w0yy,w0yz,w0zz
457 
458  real delta22,delta23,delta42,delta43
459 
460 #If #ORDER == "2"
461  declareDifferenceOrder2(u,RX)
462 #End
463 #If #ORDER == "4"
464  declareDifferenceOrder4(u,RX)
465 #End
466  ! --- begin statement functions
467  rx(i1,i2,i3)=rsxy(i1,i2,i3,0,0)
468  ry(i1,i2,i3)=rsxy(i1,i2,i3,0,1)
469  rz(i1,i2,i3)=rsxy(i1,i2,i3,0,2)
470  sx(i1,i2,i3)=rsxy(i1,i2,i3,1,0)
471  sy(i1,i2,i3)=rsxy(i1,i2,i3,1,1)
472  sz(i1,i2,i3)=rsxy(i1,i2,i3,1,2)
473  tx(i1,i2,i3)=rsxy(i1,i2,i3,2,0)
474  ty(i1,i2,i3)=rsxy(i1,i2,i3,2,1)
475  tz(i1,i2,i3)=rsxy(i1,i2,i3,2,2)
476 
477 c The next macro call will define the difference approximation statement functions
478 #If #ORDER == "2"
479  defineDifferenceOrder2Components1(u,RX)
480 #End
481 #If #ORDER == "4"
482  defineDifferenceOrder4Components1(u,RX)
483 #End
484 
485 c --- For 2nd order 2D artificial diffusion ---
486  delta22(c)= (u(i1+1,i2,i3,c)-4.*u(i1,i2,i3,c)+u(i1-1,i2,i3,c) \
487  +u(i1,i2+1,i3,c) +u(i1,i2-1,i3,c))
488 c --- For 2nd order 3D artificial diffusion ---
489  delta23(c)= \
490  (u(i1+1,i2,i3,c)-6.*u(i1,i2,i3,c)+u(i1-1,i2,i3,c) \
491  +u(i1,i2+1,i3,c) +u(i1,i2-1,i3,c) \
492  +u(i1,i2,i3+1,c) +u(i1,i2,i3-1,c))
493 c ---For fourth-order artificial diffusion in 2D
494  delta42(c)= \
495  ( -u(i1+2,i2,i3,c)-u(i1-2,i2,i3,c) \
496  -u(i1,i2+2,i3,c)-u(i1,i2-2,i3,c) \
497  +4.*(u(i1+1,i2,i3,c)+u(i1-1,i2,i3,c) \
498  +u(i1,i2+1,i3,c)+u(i1,i2-1,i3,c)) \
499  -12.*u(i1,i2,i3,c) )
500 c ---For fourth-order artificial diffusion in 3D
501  delta43(c)= \
502  ( -u(i1+2,i2,i3,c)-u(i1-2,i2,i3,c) \
503  -u(i1,i2+2,i3,c)-u(i1,i2-2,i3,c) \
504  -u(i1,i2,i3+2,c)-u(i1,i2,i3-2,c) \
505  +4.*(u(i1+1,i2,i3,c)+u(i1-1,i2,i3,c) \
506  +u(i1,i2+1,i3,c)+u(i1,i2-1,i3,c) \
507  +u(i1,i2,i3+1,c)+u(i1,i2,i3-1,c)) \
508  -18.*u(i1,i2,i3,c) )
509 
510 c --- end statement functions
511 
512  ierr=0
513  ! write(*,'("Inside insdt: gridType=",i2)') gridType
514 
515  pc =ipar(0)
516  uc =ipar(1)
517  vc =ipar(2)
518  wc =ipar(3)
519  nc =ipar(4)
520  sc =ipar(5)
521  tc =ipar(6) ! **new**
522  grid =ipar(7)
523  orderOfAccuracy =ipar(8)
524  gridIsMoving =ipar(9)
525  useWhereMask =ipar(10)
526  gridIsImplicit =ipar(11)
527  implicitMethod =ipar(12)
528  implicitOption =ipar(13)
529  isAxisymmetric =ipar(14)
530  use2ndOrderAD =ipar(15)
531  use4thOrderAD =ipar(16)
532  advectPassiveScalar=ipar(17)
533  gridType =ipar(18)
534  turbulenceModel =ipar(19)
535  pdeModel =ipar(20)
536  vsc =ipar(21)
537  rc =ipar(22)
538 
539  dr(0) =rpar(0)
540  dr(1) =rpar(1)
541  dr(2) =rpar(2)
542  dx(0) =rpar(3)
543  dx(1) =rpar(4)
544  dx(2) =rpar(5)
545  nu =rpar(6)
546  ad21 =rpar(7)
547  ad22 =rpar(8)
548  ad41 =rpar(9)
549  ad42 =rpar(10)
550  nuPassiveScalar =rpar(11)
551  adcPassiveScalar =rpar(12)
552  ad21n =rpar(13)
553  ad22n =rpar(14)
554  ad41n =rpar(15)
555  ad42n =rpar(16)
556 
557  ! nuVP =rpar(24) ! for visco-plastic
558  ! etaVP =rpar(25)
559  ! yieldStressVP =rpar(26)
560  ! exponentVP =rpar(27)
561  ! epsVP =rpar(28)
562 
563 ! write(*,'(" insdt: eta,yield,exp,eps=",4e10.2)') etaVP,yieldStressVP,exponentVP,epsVP
564 
565  kc=nc
566  ec=kc+1
567 
568  if( orderOfAccuracy.ne.2 .and. orderOfAccuracy.ne.4 )then
569  write(*,'("insdt:ERROR orderOfAccuracy=",i6)') orderOfAccuracy
570  stop 1
571  end if
572  if( gridType.ne.rectangular .and. gridType.ne.curvilinear )then
573  write(*,'("insdt:ERROR gridType=",i6)') gridType
574  stop 2
575  end if
576  if( uc.lt.0 .or. vc.lt.0 .or. (nd.eq.3 .and. wc.lt.0) )then
577  write(*,'("insdt:ERROR uc,vc,ws=",3i6)') uc,vc,wc
578  stop 4
579  end if
580 
581 c write(*,'("insdt: turbulenceModel=",2i6)') turbulenceModel
582 c write(*,'("insdt: nd,uc,vc,wc,kc=",2i6)') nd,uc,vc,wc,kc
583 
584  if( turbulenceModel.eq.kEpsilon .and. (kc.lt.uc+nd .or. kc.gt.1000) )then
585  write(*,'("insdt:ERROR in kc: nd,uc,vc,wc,kc=",2i6)') nd,uc,vc,wc,kc
586  stop 5
587  end if
588 
589 
590 c ** these are needed by self-adjoint terms **fix**
591  dxi=1./dx(0)
592  dyi=1./dx(1)
593  dzi=1./dx(2)
594  dri=1./dr(0)
595  dsi=1./dr(1)
596  dti=1./dr(2)
597  dr2i=1./(2.*dr(0))
598  ds2i=1./(2.*dr(1))
599  dt2i=1./(2.*dr(2))
600 
601  if( turbulenceModel.eq.spalartAllmaras )then
602  call getSpalartAllmarasParameters(cb1, cb2, cv1, sigma, sigmai, kappa, cw1, cw2, cw3, cw3e6, \
603  cv1e3, cd0, cr0)
604  else if( turbulenceModel.eq.kEpsilon )then
605 
606  ! write(*,'(" insdt: k-epsilon: nc,kc,ec=",3i3)') nc,kc,ec
607 
608  call getKEpsilonParameters( cMu,cEps1,cEps2,sigmaEpsI,sigmaKI )
609  ! write(*,'(" insdt: cMu,cEps1,cEps2,sigmaEpsI,sigmaKI=",5f8.3)') cMu,cEps1,cEps2,sigmaEpsI,sigmaKI
610 
611  else if( turbulenceModel.eq.largeEddySimulation )then
612  ! do nothing
613  else if( turbulenceModel.ne.noTurbulenceModel )then
614  stop 88
615  end if
616 
617  adc=adcPassiveScalar ! coefficient of linear artificial diffusion
618  cd22=ad22/(nd**2)
619  cd42=ad42/(nd**2)
620 
621 c *********************************
622 c ********MAIN LOOPS***************
623 c *********************************
624  assignEquations(DIM,ORDER)
625 
626  return
627  end
628 #endMacro
629 
630 c
631 c : empty version for linking when we do not want an option
632 c
633 #beginMacro INSDT_NULL(NAME,DIM,ORDER)
634  subroutine NAME(nd,n1a,n1b,n2a,n2b,n3a,n3b,nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b,\
635  mask,xy,rsxy,radiusInverse, u,uu, ut,uti,gv,dw, bc, ipar, rpar, ierr )
636 c======================================================================
637 c EMPTY VERSION for Linking without this Capability
638 c
639 c Compute du/dt for the incompressible NS on rectangular grids
640 c OPTIMIZED version for rectangular grids.
641 c nd : number of space dimensions
642 c
643 c gv : gridVelocity for moving grids
644 c uu : for moving grids uu is a workspace to hold u-gv, otherwise uu==u
645 c dw : distance to the wall for some turbulence models
646 c======================================================================
647  implicit none
648  integer nd, n1a,n1b,n2a,n2b,n3a,n3b,nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b
649 
650  real u(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
651  real uu(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
652  real ut(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
653  real uti(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
654  real gv(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
655  real dw(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
656  real xy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
657  real rsxy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1,0:nd-1)
658  real radiusInverse(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
659 
660  integer mask(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
661  integer bc(0:1,0:2),ierr
662 
663  integer ipar(0:*)
664  real rpar(0:*)
665 
666  write(*,'("ERROR: NULL version of subroutine NAME called")')
667  write(*,'(" You may have to turn on an option in the Makefile.")')
668 
669  stop 1080
670 
671  return
672  end
673 #endMacro
674 
675 
676 #beginMacro buildFile(NAME,DIM,ORDER)
677 #beginFile src/NAME.f
678  INSDT(NAME,DIM,ORDER)
679 #endFile
680 #beginFile src/NAME ## Null.f
681  INSDT_NULL(NAME,DIM,ORDER)
682 #endFile
683 #endMacro
684 
685 
Generated on Fri Jan 4 2013 10:19:21 for CG by   doxygen 1.8.3