CG  Version 25
lineSolveINS.h
Go to the documentation of this file.
1 !
2 ! --- Macros to define INS and Temperature line-solve functions:---
3 ! fillEquationsRectangularGridINS
4 ! fillEquationsCurvilinearGridINS
5 ! fillEquationsRectangularGridTemperature
6 ! fillEquationsCurvilinearGridTemperature
7 ! computeResidualINS
8 
9 ! This file is included in insLineSolveNew.bf
10 
11 
12 c ===================================================================================
13 c INS: Incompressible Navier Stokes
14 c dir: 0,1,2
15 c====================================================================================
16 #beginMacro fillEquationsRectangularGridINS(dir)
17 
18  ! ****** Incompressible NS, No turbulence model *****
19  if( use4thOrderAD.eq.1 )then
20  write(*,*) 'insLineSolve: 4th order diss not finished'
21  stop 7654
22  end if
23 
24  ! set default values for no 2nd order artificial diffusion:
25  cdm=0.
26  cdDiag=0.
27  cdp=0.
28 
29 ! INS - RHS forcing for rectangular grids, directions=0,1,2 (do NOT include grad(p) terms,
30 ! since then macro is valid for m=uc,vc,wc)
31 #If #dir == "0"
32  #defineMacro fr2d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(vc)*uy2(m)+nu*uyy0(m) \
33  -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
34  #defineMacro fr3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(vc)*uy2(m)-uu(wc)*uz2(m)+nu*(uyy0(m)+uzz0(m))\
35  -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
36 #Elif #dir == "1"
37  #defineMacro fr2d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)+nu*uxx0(m) \
38  -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
39  #defineMacro fr3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)-uu(wc)*uz2(m)+nu*(uxx0(m)+uzz0(m)) -\
40  thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
41 #Else
42  #defineMacro fr2d(m) (0.)
43  #defineMacro fr3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)-uu(vc)*uy2(m)+nu*(uxx0(m)+uyy0(m)) \
44  -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
45 #End
46 
47 if( nd.eq.2 )then
48 
49  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
50  defineDerivativeMacros(2,2,rectangular)
51  beginLoops()
52  if( mask(i1,i2,i3).gt.0 )then
53  if( use2ndOrderAD.eq.1 )then
54  defineArtificialDiffusionCoefficients(2,dir,R,)
55  end if
56  if( computeMatrix.eq.1 )then
57  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir)-nu*dxvsqi(dir) -cdm
58  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*nu*(dxvsqi(0)+dxvsqi(1)) +cdDiag
59  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir)-nu*dxvsqi(dir) -cdp
60  end if
61  if( computeRHS.eq.1 )then
62  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+fr2d(uc)-ux2(pc)
63  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+fr2d(vc)-uy2(pc)
64  if( use2ndOrderAD.eq.1 )then
65  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+ adE ## dir(i1,i2,i3,uc)
66  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+ adE ## dir(i1,i2,i3,vc)
67  end if
68  end if
69  else
70  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
71  am(i1,i2,i3)=0.
72  bm(i1,i2,i3)=1.
73  cm(i1,i2,i3)=0.
74  end if
75  if( computeRHS.eq.1 )then
76  f(i1,i2,i3,fcu)=uu(uc)
77  f(i1,i2,i3,fcv)=uu(vc)
78  end if
79  end if
80  endLoops()
81 
82 else if( nd.eq.3 )then
83 
84  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE)
85  defineDerivativeMacros(3,2,rectangular)
86 
87  beginLoops()
88  if( mask(i1,i2,i3).gt.0 )then
89  if( use2ndOrderAD.eq.1 )then
90  defineArtificialDiffusionCoefficients(3,dir,R,)
91  end if
92  if( computeMatrix.eq.1 )then
93  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir)-nu*dxvsqi(dir) -cdm
94  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*nu*(dxvsqi(0)+dxvsqi(1)+dxvsqi(2)) +cdDiag
95  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir)-nu*dxvsqi(dir) -cdp
96  end if
97  if( computeRHS.eq.1 )then
98  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+fr3d(uc)-ux2(pc)
99  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+fr3d(vc)-uy2(pc)
100  f(i1,i2,i3,fcw)=f(i1,i2,i3,fcw)+fr3d(wc)-uz2(pc)
101  if( use2ndOrderAD.eq.1 )then
102  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+adE3d ## dir(i1,i2,i3,uc)
103  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+adE3d ## dir(i1,i2,i3,vc)
104  f(i1,i2,i3,fcw)=f(i1,i2,i3,fcw)+adE3d ## dir(i1,i2,i3,wc)
105  end if
106  end if
107  else
108  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
109  am(i1,i2,i3)=0.
110  bm(i1,i2,i3)=1.
111  cm(i1,i2,i3)=0.
112  end if
113  if( computeRHS.eq.1 )then
114  f(i1,i2,i3,fcu)=uu(uc)
115  f(i1,i2,i3,fcv)=uu(vc)
116  f(i1,i2,i3,fcw)=uu(wc)
117  end if
118  end if
119  endLoops()
120 
121 else
122  stop 888 ! unexpected value for nd
123 end if
124 
125 #endMacro
126 
127 
128 c ===================================================================================
129 c ****** Incompressible NS, No turbulence model *****
130 c dir: 0,1,2
131 c====================================================================================
132 #beginMacro fillEquationsCurvilinearGridINS(dir)
133 
134  if( use4thOrderAD.eq.1 )then
135  write(*,*) 'insLineSolve: 4th order diss not finished'
136  stop 7655
137  end if
138 
139  dirp1=mod(dir+1,nd)
140  dirp2=mod(dir+2,nd)
141 
142  ! INS - RHS forcing for curvilinear grids, directions=0,1,2 (do NOT include grad(p) terms)
143 #If #dir == "0"
144  #defineMacro fc2d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
145  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)+nu*(rxx(1,0)+ryy(1,1)))*us(m)+\
146  nu*((rxi(1,0)**2+rxi(1,1)**2)*uss0(m)+\
147  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1))*urs(m)) -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
148  #defineMacro fc3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
149  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)-uu(wc)*rxi(1,2)+nu*(rxx3(1,0)+rxy3(1,1)+rxz3(1,2)))*us(m)+\
150  (-uu(uc)*rxi(2,0)-uu(vc)*rxi(2,1)-uu(wc)*rxi(2,2)+nu*(rxx3(2,0)+rxy3(2,1)+rxz3(2,2)))*ut(m)+\
151  nu*( (rxi(1,0)**2+rxi(1,1)**2+rxi(1,2)**2)*uss0(m)+\
152  (rxi(2,0)**2+rxi(2,1)**2+rxi(2,2)**2)*utt0(m)+\
153  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
154  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
155  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
156  ) -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
157 #Elif #dir == "1"
158  #defineMacro fc2d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
159  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)+nu*(rxx(0,0)+ryy(0,1)))*ur(m)+\
160  nu*((rxi(0,0)**2+rxi(0,1)**2)*urr0(m)+\
161  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1))*urs(m)) -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
162  #defineMacro fc3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
163  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)-uu(wc)*rxi(0,2)+nu*(rxx3(0,0)+rxy3(0,1)+rxz3(0,2)))*ur(m)+\
164  (-uu(uc)*rxi(2,0)-uu(vc)*rxi(2,1)-uu(wc)*rxi(2,2)+nu*(rxx3(2,0)+rxy3(2,1)+rxz3(2,2)))*ut(m)+\
165  nu*( (rxi(0,0)**2+rxi(0,1)**2+rxi(0,2)**2)*urr0(m)+\
166  (rxi(2,0)**2+rxi(2,1)**2+rxi(2,2)**2)*utt0(m)+\
167  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
168  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
169  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
170  ) -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
171 #Else
172  #defineMacro fc2d(m) (0.)
173  #defineMacro fc3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
174  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)-uu(wc)*rxi(0,2)+nu*(rxx3(0,0)+rxy3(0,1)+rxz3(0,2)))*ur(m)+\
175  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)-uu(wc)*rxi(1,2)+nu*(rxx3(1,0)+rxy3(1,1)+rxz3(1,2)))*us(m)+\
176  nu*( (rxi(0,0)**2+rxi(0,1)**2+rxi(0,2)**2)*urr0(m)+\
177  (rxi(1,0)**2+rxi(1,1)**2+rxi(1,2)**2)*uss0(m)+\
178  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
179  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
180  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
181  ) -thermalExpansivity*gravity(m-uc)*u(i1,i2,i3,tc))
182 #End
183 
184  ! set default values for no 2nd order artificial diffusion:
185  cdm=0.
186  cdDiag=0.
187  cdp=0.
188 
189 if( nd.eq.2 )then
190  defineDerivativeMacros(2,2,curvilinear)
191 
192  beginLoops()
193  if( mask(i1,i2,i3).gt.0 )then
194  if( use2ndOrderAD.eq.1 )then
195  defineArtificialDiffusionCoefficients(2,dir,C,)
196  end if
197  if( computeMatrix.eq.1 )then
198  t1=(uu(uc)*rxi(dir,0)+uu(vc)*rxi(dir,1)-nu*(rxx(dir,0)+rxy(dir,1)))*drv2i(dir)
199  t2=nu*(rxi(dir,0)**2+rxi(dir,1)**2)*drvsqi(dir)
200  am(i1,i2,i3)= -t1-t2 -cdm
201  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*(t2+nu*(rxi(dirp1,0)**2+rxi(dirp1,1)**2)*drvsqi(dirp1) )+cdDiag
202  cm(i1,i2,i3)= t1-t2 -cdp
203  end if
204  if( computeRHS.eq.1 )then
205  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+fc2d(uc)-ux2c(pc)
206  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+fc2d(vc)-uy2c(pc)
207  if( use2ndOrderAD.eq.1 )then
208  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+adE ## dir(i1,i2,i3,uc)
209  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+adE ## dir(i1,i2,i3,vc)
210  end if
211  end if
212  else ! for interpolation points or unused:
213  if( computeMatrix.eq.1 )then
214  am(i1,i2,i3)=0.
215  bm(i1,i2,i3)=1.
216  cm(i1,i2,i3)=0.
217  end if
218  if( computeRHS.eq.1 )then
219  f(i1,i2,i3,fcu)=uu(uc)
220  f(i1,i2,i3,fcv)=uu(vc)
221  end if
222  end if
223  endLoops()
224 
225 else if( nd.eq.3 )then
226 
227  defineDerivativeMacros(3,2,curvilinear)
228 
229  beginLoops()
230  if( mask(i1,i2,i3).gt.0 )then
231  if( use2ndOrderAD.eq.1 )then
232  defineArtificialDiffusionCoefficients(3,dir,C,)
233  end if
234  if( computeMatrix.eq.1 )then
235  t1=(uu(uc)*rxi(dir,0)+uu(vc)*rxi(dir,1)+uu(wc)*rxi(dir,2)-nu*(rxx3(dir,0)+rxy3(dir,1)+rxz3(dir,2)))*drv2i(dir)
236  t2=nu*(rxi(dir,0)**2+rxi(dir,1)**2+rxi(dir,2)**2)*drvsqi(dir)
237  am(i1,i2,i3)= -t1-t2 -cdm
238  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3)\
239  +2.*(t2+nu*( (rxi(dirp1,0)**2+rxi(dirp1,1)**2+rxi(dirp1,2)**2)*drvsqi(dirp1)+\
240  (rxi(dirp2,0)**2+rxi(dirp2,1)**2+rxi(dirp2,2)**2)*drvsqi(dirp2) )) +cdDiag
241  cm(i1,i2,i3)= t1-t2 -cdp
242  end if
243  if( computeRHS.eq.1 )then
244  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+fc3d(uc)-ux3c(pc)
245  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+fc3d(vc)-uy3c(pc)
246  f(i1,i2,i3,fcw)=f(i1,i2,i3,fcw)+fc3d(wc)-uz3c(pc)
247  if( use2ndOrderAD.eq.1 )then
248  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+adE3d ## dir(i1,i2,i3,uc)
249  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+adE3d ## dir(i1,i2,i3,vc)
250  f(i1,i2,i3,fcw)=f(i1,i2,i3,fcw)+adE3d ## dir(i1,i2,i3,wc)
251  end if
252  end if
253  else ! for interpolation points or unused:
254  if( computeMatrix.eq.1 )then
255  am(i1,i2,i3)=0.
256  bm(i1,i2,i3)=1.
257  cm(i1,i2,i3)=0.
258  end if
259  if( computeRHS.eq.1 )then
260  f(i1,i2,i3,fcu)=uu(uc)
261  f(i1,i2,i3,fcv)=uu(vc)
262  f(i1,i2,i3,fcw)=uu(wc)
263  end if
264  end if
265  endLoops()
266 
267 else
268  stop 222 ! unexpected value for nd
269 end if
270 
271 
272 #endMacro
273 
274 
275 c ===================================================================================
276 c INS: Incompressible Navier Stokes Temperature Equation
277 c dir: 0,1,2
278 c====================================================================================
279 #beginMacro fillEquationsRectangularGridTemperature(dir)
280 
281  ! write(*,*) 'new: fillEquationsRectangularGridTemperature'
282 
283  ! ****** Temperature Equation for INS *****
284  if( use4thOrderAD.eq.1 )then
285  write(*,*) 'insLineSolve: T : 4th order diss not finished'
286  stop 7654
287  end if
288 
289  ! set default values for no 2nd order artificial diffusion:
290  cdm=0.
291  cdDiag=0.
292  cdp=0.
293 
294 #If #dir == "0"
295  #defineMacro frt2d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(vc)*uy2(m)+kThermal*uyy0(m))
296  #defineMacro frt3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(vc)*uy2(m)-uu(wc)*uz2(m)+kThermal*(uyy0(m)+uzz0(m)))
297 #Elif #dir == "1"
298  #defineMacro frt2d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)+kThermal*uxx0(m))
299  #defineMacro frt3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)-uu(wc)*uz2(m)+kThermal*(uxx0(m)+uzz0(m)))
300 #Else
301  #defineMacro frt2d(m) (0.)
302  #defineMacro frt3d(m) (uu(m)*dtScale/dt(i1,i2,i3) -uu(uc)*ux2(m)-uu(vc)*uy2(m)+kThermal*(uxx0(m)+uyy0(m)))
303 #End
304 
305 if( nd.eq.2 )then
306 
307  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
308  defineDerivativeMacros(2,2,rectangular)
309  beginLoops()
310  if( mask(i1,i2,i3).gt.0 )then
311  if( use2ndOrderAD.eq.1 )then
312  defineArtificialDiffusionCoefficients(2,dir,R,)
313  end if
314  if( computeMatrix.eq.1 )then
315  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir)-kThermal*dxvsqi(dir) -cdm
316  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*kThermal*(dxvsqi(0)+dxvsqi(1)) +cdDiag
317  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir)-kThermal*dxvsqi(dir) -cdp
318  end if
319  if( computeRHS.eq.1 )then
320  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+frt2d(tc)
321  if( use2ndOrderAD.eq.1 )then
322  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+ adE ## dir(i1,i2,i3,tc)
323  end if
324  end if
325  else
326  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
327  am(i1,i2,i3)=0.
328  bm(i1,i2,i3)=1.
329  cm(i1,i2,i3)=0.
330  end if
331  if( computeRHS.eq.1 )then
332  f(i1,i2,i3,fct)=uu(tc)
333  end if
334  end if
335  endLoops()
336 
337 else if( nd.eq.3 )then
338 
339  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE)
340  defineDerivativeMacros(3,2,rectangular)
341 
342  beginLoops()
343  if( mask(i1,i2,i3).gt.0 )then
344  if( use2ndOrderAD.eq.1 )then
345  defineArtificialDiffusionCoefficients(3,dir,R,)
346  end if
347  if( computeMatrix.eq.1 )then
348  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir)-kThermal*dxvsqi(dir) -cdm
349  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*kThermal*(dxvsqi(0)+dxvsqi(1)+dxvsqi(2)) +cdDiag
350  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir)-kThermal*dxvsqi(dir) -cdp
351  end if
352  if( computeRHS.eq.1 )then
353  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+frt3d(tc)
354  if( use2ndOrderAD.eq.1 )then
355  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+adE3d ## dir(i1,i2,i3,tc)
356  end if
357  end if
358  else
359  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
360  am(i1,i2,i3)=0.
361  bm(i1,i2,i3)=1.
362  cm(i1,i2,i3)=0.
363  end if
364  if( computeRHS.eq.1 )then
365  f(i1,i2,i3,fct)=uu(tc)
366  end if
367  end if
368  endLoops()
369 
370 else
371  stop 888 ! unexpected value for nd
372 end if
373 
374 #endMacro
375 
376 
377 c ===================================================================================
378 c ****** Temperature Equation for INS *****
379 c dir: 0,1,2
380 c====================================================================================
381 #beginMacro fillEquationsCurvilinearGridTemperature(dir)
382 
383  ! write(*,*) 'new: fillEquationsCurvilinearGridTemperature'
384 
385  if( use4thOrderAD.eq.1 )then
386  write(*,*) 'insLineSolve: T : 4th order diss not finished'
387  stop 7655
388  end if
389 
390  dirp1=mod(dir+1,nd)
391  dirp2=mod(dir+2,nd)
392 
393 #If #dir == "0"
394  #defineMacro fct2d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
395  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)+kThermal*(rxx(1,0)+ryy(1,1)))*us(m)+\
396  kThermal*((rxi(1,0)**2+rxi(1,1)**2)*uss0(m)+\
397  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1))*urs(m)) )
398  #defineMacro fct3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
399  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)-uu(wc)*rxi(1,2)+kThermal*(rxx3(1,0)+rxy3(1,1)+rxz3(1,2)))*us(m)+\
400  (-uu(uc)*rxi(2,0)-uu(vc)*rxi(2,1)-uu(wc)*rxi(2,2)+kThermal*(rxx3(2,0)+rxy3(2,1)+rxz3(2,2)))*ut(m)+\
401  kThermal*( (rxi(1,0)**2+rxi(1,1)**2+rxi(1,2)**2)*uss0(m)+\
402  (rxi(2,0)**2+rxi(2,1)**2+rxi(2,2)**2)*utt0(m)+\
403  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
404  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
405  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
406  ) )
407 #Elif #dir == "1"
408  #defineMacro fct2d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
409  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)+kThermal*(rxx(0,0)+ryy(0,1)))*ur(m)+\
410  kThermal*((rxi(0,0)**2+rxi(0,1)**2)*urr0(m)+\
411  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1))*urs(m)) )
412  #defineMacro fct3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
413  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)-uu(wc)*rxi(0,2)+kThermal*(rxx3(0,0)+rxy3(0,1)+rxz3(0,2)))*ur(m)+\
414  (-uu(uc)*rxi(2,0)-uu(vc)*rxi(2,1)-uu(wc)*rxi(2,2)+kThermal*(rxx3(2,0)+rxy3(2,1)+rxz3(2,2)))*ut(m)+\
415  kThermal*( (rxi(0,0)**2+rxi(0,1)**2+rxi(0,2)**2)*urr0(m)+\
416  (rxi(2,0)**2+rxi(2,1)**2+rxi(2,2)**2)*utt0(m)+\
417  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
418  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
419  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
420  ) )
421 #Else
422  #defineMacro fct2d(m) (0.)
423  #defineMacro fct3d(m) (uu(m)*dtScale/dt(i1,i2,i3) +\
424  (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1)-uu(wc)*rxi(0,2)+kThermal*(rxx3(0,0)+rxy3(0,1)+rxz3(0,2)))*ur(m)+\
425  (-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1)-uu(wc)*rxi(1,2)+kThermal*(rxx3(1,0)+rxy3(1,1)+rxz3(1,2)))*us(m)+\
426  kThermal*( (rxi(0,0)**2+rxi(0,1)**2+rxi(0,2)**2)*urr0(m)+\
427  (rxi(1,0)**2+rxi(1,1)**2+rxi(1,2)**2)*uss0(m)+\
428  2.*(rxi(0,0)*rxi(1,0)+rxi(0,1)*rxi(1,1)+rxi(0,2)*rxi(1,2))*urs(m)+\
429  2.*(rxi(0,0)*rxi(2,0)+rxi(0,1)*rxi(2,1)+rxi(0,2)*rxi(2,2))*urt(m)+\
430  2.*(rxi(1,0)*rxi(2,0)+rxi(1,1)*rxi(2,1)+rxi(1,2)*rxi(2,2))*ust(m)\
431  ) )
432 #End
433 
434  ! set default values for no 2nd order artificial diffusion:
435  cdm=0.
436  cdDiag=0.
437  cdp=0.
438 
439 if( nd.eq.2 )then
440  defineDerivativeMacros(2,2,curvilinear)
441 
442  beginLoops()
443  if( mask(i1,i2,i3).gt.0 )then
444  if( use2ndOrderAD.eq.1 )then
445  defineArtificialDiffusionCoefficients(2,dir,C,)
446  end if
447  if( computeMatrix.eq.1 )then
448  t1=(uu(uc)*rxi(dir,0)+uu(vc)*rxi(dir,1)-kThermal*(rxx(dir,0)+rxy(dir,1)))*drv2i(dir)
449  t2=kThermal*(rxi(dir,0)**2+rxi(dir,1)**2)*drvsqi(dir)
450  am(i1,i2,i3)= -t1-t2 -cdm
451  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +2.*(t2+kThermal*(rxi(dirp1,0)**2+rxi(dirp1,1)**2)*drvsqi(dirp1) )+cdDiag
452  cm(i1,i2,i3)= t1-t2 -cdp
453  end if
454  if( computeRHS.eq.1 )then
455  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+fct2d(tc)
456  if( use2ndOrderAD.eq.1 )then
457  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+adE ## dir(i1,i2,i3,tc)
458  end if
459  end if
460  else ! for interpolation points or unused:
461  if( computeMatrix.eq.1 )then
462  am(i1,i2,i3)=0.
463  bm(i1,i2,i3)=1.
464  cm(i1,i2,i3)=0.
465  end if
466  if( computeRHS.eq.1 )then
467  f(i1,i2,i3,fct)=uu(tc)
468  end if
469  end if
470  endLoops()
471 
472 else if( nd.eq.3 )then
473 
474  defineDerivativeMacros(3,2,curvilinear)
475 
476  beginLoops()
477  if( mask(i1,i2,i3).gt.0 )then
478  if( use2ndOrderAD.eq.1 )then
479  defineArtificialDiffusionCoefficients(3,dir,C,)
480  end if
481  if( computeMatrix.eq.1 )then
482  t1=(uu(uc)*rxi(dir,0)+uu(vc)*rxi(dir,1)+uu(wc)*rxi(dir,2)-kThermal*(rxx3(dir,0)+rxy3(dir,1)+rxz3(dir,2)))*drv2i(dir)
483  t2=kThermal*(rxi(dir,0)**2+rxi(dir,1)**2+rxi(dir,2)**2)*drvsqi(dir)
484  am(i1,i2,i3)= -t1-t2 -cdm
485  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3)\
486  +2.*(t2+kThermal*( (rxi(dirp1,0)**2+rxi(dirp1,1)**2+rxi(dirp1,2)**2)*drvsqi(dirp1)+\
487  (rxi(dirp2,0)**2+rxi(dirp2,1)**2+rxi(dirp2,2)**2)*drvsqi(dirp2) )) +cdDiag
488  cm(i1,i2,i3)= t1-t2 -cdp
489  end if
490  if( computeRHS.eq.1 )then
491  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+fct3d(tc)
492  if( use2ndOrderAD.eq.1 )then
493  f(i1,i2,i3,fct)=f(i1,i2,i3,fct)+adE3d ## dir(i1,i2,i3,tc)
494  end if
495  end if
496  else ! for interpolation points or unused:
497  if( computeMatrix.eq.1 )then
498  am(i1,i2,i3)=0.
499  bm(i1,i2,i3)=1.
500  cm(i1,i2,i3)=0.
501  end if
502  if( computeRHS.eq.1 )then
503  f(i1,i2,i3,fct)=uu(tc)
504  end if
505  end if
506  endLoops()
507 
508 else
509  stop 222 ! unexpected value for nd
510 end if
511 
512 
513 #endMacro
514 
515 
516 
517 c==================================================================================
518 c Residual Computation for INS: incompressible Navier Stokes (including Boussinesq)
519 c
520 c Macro args:
521 c GRIDTYPE: rectangular, curvilinear
522 c====================================================================================
523 #beginMacro computeResidualINS(GRIDTYPE)
524 
525  ! write(*,*) 'new: computeResidualINS'
526 
527  if( nd.eq.2 )then
528 
529  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
530  ! * defineDerivativeMacros(DIM,2,GRIDTYPE)
531 
532 
533  beginLoops()
534  if( mask(i1,i2,i3).gt.0 )then
535 
536  #If #GRIDTYPE == "rectangular"
537  ! ************ INS 2d rectangular case ********************
538 
539  ! u.t + u.grad(u) + p.x = nu Delta(u)
540  ! v.t + u.grad(v) + p.y = nu Delta(v)
541  residual(i1,i2,i3,uc)=f(i1,i2,i3,uc)-u(i1,i2,i3,uc)*ux2(uc)-u(i1,i2,i3,vc)*uy2(uc)-ux2(pc)+nu*lap2d2(uc)\
542  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc)
543 
544  residual(i1,i2,i3,vc)=f(i1,i2,i3,vc)-u(i1,i2,i3,uc)*ux2(vc)-u(i1,i2,i3,vc)*uy2(vc)-uy2(pc)+nu*lap2d2(vc)\
545  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc)
546 
547  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
548  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint2dR(i1,i2,i3,uc)
549  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint2dR(i1,i2,i3,vc)
550  else if( use4thOrderAD.eq.1 )then
551  ! compute adc2, adc4:
552  getDerivativesAndDissipation(2,2,GRIDTYPE)
553  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+ad2(adc2,uc)+ad4(adc4,uc)
554  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+ad2(adc2,vc)+ad4(adc4,vc)
555  end if
556 
557  if( computeTemperature.ne.0 )then
558  residual(i1,i2,i3,tc)=f(i1,i2,i3,tc)-u(i1,i2,i3,uc)*ux2(tc)-u(i1,i2,i3,vc)*uy2(tc)+kThermal*lap2d2(tc)
559  ! --- artificial dissipation for T: do this for now:
560  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
561  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint2dR(i1,i2,i3,tc)
562  else if( use4thOrderAD.eq.1 )then
563  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad2(adc2,tc)+ad4(adc4,tc)
564  end if
565  end if
566 
567  #Else
568  ! ************ INS 2d curvilinear case ********************
569 
570  residual(i1,i2,i3,uc)=f(i1,i2,i3,uc)-u(i1,i2,i3,uc)*ux2c(uc)-u(i1,i2,i3,vc)*uy2c(uc)-ux2c(pc)+nu*lap2d2c(uc)\
571  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc)
572  residual(i1,i2,i3,vc)=f(i1,i2,i3,vc)-u(i1,i2,i3,uc)*ux2c(vc)-u(i1,i2,i3,vc)*uy2c(vc)-uy2c(pc)+nu*lap2d2c(vc)\
573  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc)
574 
575  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
576  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint2dC(i1,i2,i3,uc)
577  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint2dC(i1,i2,i3,vc)
578  else if( use4thOrderAD.eq.1 )then
579  ! compute adc2, adc4:
580  getDerivativesAndDissipation(2,2,GRIDTYPE)
581  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+ad2(adc2,uc)+ad4(adc4,uc)
582  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+ad2(adc2,vc)+ad4(adc4,vc)
583  end if
584 
585  if( computeTemperature.ne.0 )then
586  residual(i1,i2,i3,tc)=f(i1,i2,i3,tc)-u(i1,i2,i3,uc)*ux2c(tc)-u(i1,i2,i3,vc)*uy2c(tc)+kThermal*lap2d2c(tc)
587  ! --- artificial dissipation for T: do this for now:
588  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
589  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint2dC(i1,i2,i3,tc)
590  else if( use4thOrderAD.eq.1 )then
591  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad2(adc2,tc)+ad4(adc4,tc)
592  end if
593  end if
594 
595  #End
596 
597  end if
598  endLoops()
599 
600  else if( nd.eq.3 )then
601 
602  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
603  ! defineDerivativeMacros(DIM,2,GRIDTYPE)
604 
605 
606  beginLoops()
607  if( mask(i1,i2,i3).gt.0 )then
608 
609  #If #GRIDTYPE == "rectangular"
610 
611  ! ************ INS 3d rectangular case ********************
612 
613  residual(i1,i2,i3,uc)=f(i1,i2,i3,uc)-u(i1,i2,i3,uc)*ux2(uc)\
614  -u(i1,i2,i3,vc)*uy2(uc)\
615  -u(i1,i2,i3,wc)*uz2(uc)\
616  -ux2(pc)+nu*lap3d2(uc)\
617  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc)
618  residual(i1,i2,i3,vc)=f(i1,i2,i3,vc)-u(i1,i2,i3,uc)*ux2(vc)\
619  -u(i1,i2,i3,vc)*uy2(vc)\
620  -u(i1,i2,i3,wc)*uz2(vc)\
621  -uy2(pc)+nu*lap3d2(vc)\
622  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc)
623  residual(i1,i2,i3,wc)=f(i1,i2,i3,wc)-u(i1,i2,i3,uc)*ux2(wc)\
624  -u(i1,i2,i3,vc)*uy2(wc)\
625  -u(i1,i2,i3,wc)*uz2(wc)\
626  -uz2(pc)+nu*lap3d2(wc)\
627  -thermalExpansivity*gravity(2)*u(i1,i2,i3,tc)
628  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
629  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint3dR(i1,i2,i3,uc)
630  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint3dR(i1,i2,i3,vc)
631  residual(i1,i2,i3,wc)=residual(i1,i2,i3,wc)+adSelfAdjoint3dR(i1,i2,i3,wc)
632  else if( use4thOrderAD.eq.1 )then
633  ! compute adc2, adc4:
634  getDerivativesAndDissipation(3,2,GRIDTYPE)
635  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+ad23(adc2,uc)+ad43(adc4,uc)
636  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+ad23(adc2,vc)+ad43(adc4,vc)
637  residual(i1,i2,i3,wc)=residual(i1,i2,i3,wc)+ad23(adc2,vc)+ad43(adc4,wc)
638  end if
639 
640 
641  if( computeTemperature.ne.0 )then
642  residual(i1,i2,i3,tc)=f(i1,i2,i3,tc)-u(i1,i2,i3,uc)*ux2(tc)\
643  -u(i1,i2,i3,vc)*uy2(tc)\
644  -u(i1,i2,i3,wc)*uz2(tc)\
645  +kThermal*lap3d2(tc)
646  ! --- artificial dissipation for T: do this for now:
647  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
648  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint3dR(i1,i2,i3,tc)
649  else if( use4thOrderAD.eq.1 )then
650  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad23(adc2,tc)+ad43(adc4,tc)
651  end if
652  end if
653 
654  #Else
655  ! ************ INS 3d curvilinear case ********************
656 
657 
658  residual(i1,i2,i3,uc)=f(i1,i2,i3,uc)-u(i1,i2,i3,uc)*ux3c(uc)\
659  -u(i1,i2,i3,vc)*uy3c(uc)\
660  -u(i1,i2,i3,wc)*uz3c(uc)\
661  -ux3c(pc)+nu*lap3d2c(uc)\
662  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc)
663  residual(i1,i2,i3,vc)=f(i1,i2,i3,vc)-u(i1,i2,i3,uc)*ux3c(vc)\
664  -u(i1,i2,i3,vc)*uy3c(vc)\
665  -u(i1,i2,i3,wc)*uz3c(vc)\
666  -uy3c(pc)+nu*lap3d2c(vc)\
667  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc)
668  residual(i1,i2,i3,wc)=f(i1,i2,i3,wc)-u(i1,i2,i3,uc)*ux3c(wc)\
669  -u(i1,i2,i3,vc)*uy3c(wc)\
670  -u(i1,i2,i3,wc)*uz3c(wc)\
671  -uz3c(pc)+nu*lap3d2c(wc)\
672  -thermalExpansivity*gravity(2)*u(i1,i2,i3,tc)
673  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
674  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint3dC(i1,i2,i3,uc)
675  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint3dC(i1,i2,i3,vc)
676  residual(i1,i2,i3,wc)=residual(i1,i2,i3,wc)+adSelfAdjoint3dC(i1,i2,i3,wc)
677  else if( use4thOrderAD.eq.1 )then
678  ! compute adc2, adc4:
679  getDerivativesAndDissipation(3,2,GRIDTYPE)
680  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+ad23(adc2,uc)+ad43(adc4,uc)
681  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+ad23(adc2,vc)+ad43(adc4,vc)
682  residual(i1,i2,i3,wc)=residual(i1,i2,i3,wc)+ad23(adc2,vc)+ad43(adc4,wc)
683  end if
684 
685  if( computeTemperature.ne.0 )then
686  residual(i1,i2,i3,tc)=f(i1,i2,i3,tc)-u(i1,i2,i3,uc)*ux3c(tc)\
687  -u(i1,i2,i3,vc)*uy3c(tc)\
688  -u(i1,i2,i3,wc)*uz3c(tc)\
689  +kThermal*lap3d2c(tc)
690  ! --- artificial dissipation for T: do this for now:
691  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
692  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint3dC(i1,i2,i3,tc)
693  else if( use4thOrderAD.eq.1 )then
694  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad23(adc2,tc)+ad43(adc4,tc)
695  end if
696  end if
697 
698 
699  #End
700 
701  end if
702  endLoops()
703 
704  else
705  stop 888 ! unexpected value for nd
706  end if
707 
708 #endMacro
Generated on Fri Jan 4 2013 10:19:23 for CG by   doxygen 1.8.3