CG  Version 25
insImp.h
Go to the documentation of this file.
1 ! ****************************************************************************************************
2 ! Define macros that are used to build the full implicit matrices for the INS, VP etc. equations
3 !
4 ! This file contains generic macros that are used by the different PDEs
5 !
6 ! Used by
7 ! insImpINS.bf
8 ! insImpVP.bf
9 ! insImpBL.bf
10 ! ****************************************************************************************************
11 
12 
13 c --- See --- op/fortranCoeff/opcoeff.bf
14 c op/include/defineConservative.h
15 c --- See --- mx/src/interfaceMacros.bf <-- mixes different orders of accuracy
16 
17 c -- define bpp macros for coefficient operators (from op/src/stencilCoeff.maple)
18 #Include opStencilCoeffOrder2.h
19 #Include opStencilCoeffOrder4.h
20 #Include opStencilCoeffOrder6.h
21 ! #Include opStencilCoeffOrder8.h
22 
23 
24 c These next include file will define the macros that will define the difference approximations (in op/src)
25 c Defines getDuDx2(u,aj,ff), getDuDxx2(u,aj,ff), getDuDx3(u,aj,ff), ... etc.
26 #Include "derivMacroDefinitions.h"
27 
28 c Define
29 c defineParametricDerivativeMacros(u,dr,dx,DIM,ORDER,COMPONENTS,MAXDERIV)
30 c defines -> ur2, us2, ux2, uy2, ... (2D)
31 c ur3, us3, ut3, ux3, uy3, uz3, ... (3D)
32 #Include "defineParametricDerivMacros.h"
33 
34 ! 2D, order=6, components=1
35 ! defineParametricDerivativeMacros(u,dr,dx,DIM,ORDER,COMPONENTS,MAXDERIV)
36 
37 ! defineParametricDerivativeMacros(u,dr,dx,2,2,1,2)
38  defineParametricDerivativeMacros(rsxy,dr,dx,3,2,2,2)
39  defineParametricDerivativeMacros(rsxy,dr,dx,3,4,2,2)
40  defineParametricDerivativeMacros(rsxy,dr,dx,3,6,2,2)
41 
42  defineParametricDerivativeMacros(u,dr,dx,3,2,1,2)
43  defineParametricDerivativeMacros(ul,dr,dx,3,2,1,2)
44 
45 
46 ! Example to define orders 2,4,6:
47 ! defineParametricDerivativeMacros(u1,dr1,dx1,2,2,1,6)
48 ! defineParametricDerivativeMacros(u1,dr1,dx1,2,4,1,4)
49 ! defineParametricDerivativeMacros(u1,dr1,dx1,2,6,1,2)
50 
51 
52 ! define macros for conservative operators: (in op/src)
53 ! defines getConservativeCoeff( OPERATOR,s,coeff ), OPERATOR=divScalarGrad, ...
54 #Include "conservativeCoefficientMacros.h"
55 
56 c -- From opcoeff.bf
57 
58 #beginMacro beginLoops()
59  do i3=n3a,n3b
60  do i2=n2a,n2b
61  do i1=n1a,n1b
62  if( mask(i1,i2,i3).gt.0 )then
63 #endMacro
64 
65 #beginMacro endLoops()
66  end if
67  end do
68  end do
69  end do
70 #endMacro
71 
72 #beginMacro beginLoopsNoMask()
73  do i3=n3a,n3b
74  do i2=n2a,n2b
75  do i1=n1a,n1b
76 #endMacro
77 
78 #beginMacro endLoopsNoMask()
79  end do
80  end do
81  end do
82 #endMacro
83 
84 ! This loop will check for mask>0 and mask != interiorBoundaryPoint
85 #beginMacro beginLoopsMixedBoundary()
86  do i3=n3a,n3b
87  do i2=n2a,n2b
88  do i1=n1a,n1b
89  ! if( btest(mask(i1,i2,i3),28) )then
90  ! write(*,'("+++ Point i=(",3i5,") is an interiorBoundaryPoint")') i1,i2,i3
91  ! end if
92  if( mask(i1,i2,i3).gt.0 .and. .not.btest(mask(i1,i2,i3),28) )then
93 #endMacro
94 
95 #beginMacro beginMatrixLoops(m1,m2,m3)
96  do m3=-halfWidth3,halfWidth3
97  do m2=-halfWidth,halfWidth
98  do m1=-halfWidth,halfWidth
99 #endMacro
100 
101 #beginMacro endMatrixLoops()
102  end do
103  end do
104  end do
105 #endMacro
106 
107 ! ======================================================================================================
108 ! Add the local matrix operator opCoeff to the global matrix coeff in equation "e" and component "c"
109 ! ======================================================================================================
110 #beginMacro addCoeff(c,e,coeff,opCoeff)
111 #If $DIM == 2
112  do m2=-halfWidth,halfWidth
113  do m1=-halfWidth,halfWidth
114  coeff(mce2(m1,m2,0,c,e),i1,i2,i3)=coeff(mce2(m1,m2,0,c,e),i1,i2,i3)+(opCoeff(ma2(m1,m2,0)))
115  end do
116  end do
117 #Else
118  do m3=-halfWidth,halfWidth
119  do m2=-halfWidth,halfWidth
120  do m1=-halfWidth,halfWidth
121  coeff(mce3(m1,m2,m3,c,e),i1,i2,i3)=coeff(mce3(m1,m2,m3,c,e),i1,i2,i3)+(opCoeff(ma3(m1,m2,m3)))
122  end do
123  end do
124  end do
125 #End
126 #endMacro
127 
128 
129 ! ======================================================================================================
130 ! Assign and add up to 10 different local matrix operators to the global matrix coeff
131 ! in equation "e" and component "c"
132 !
133 ! (Leave final unused arguments empty)
134 ! ======================================================================================================
135 #beginMacro setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9,op10)
136  do m3=-halfWidth3,halfWidth3
137  do m2=-halfWidth,halfWidth
138  do m1=-halfWidth,halfWidth
139 #If #op10 ne ""
140  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
141  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))+(op9(MA(m1,m2,m3)))+(op10(MA(m1,m2,m3)))
142 #Elif #op9 ne ""
143  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
144  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))+(op9(MA(m1,m2,m3)))
145 #Elif #op8 ne ""
146  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
147  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))
148 #Elif #op7 ne ""
149  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
150  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))
151 #Elif #op6 ne ""
152  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
153  (op6(MA(m1,m2,m3)))
154 #Elif #op5 ne ""
155  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))
156 #Elif #op4 ne ""
157  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))
158 #Elif #op3 ne ""
159  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))
160 #Elif #op2 ne ""
161  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))
162 #Else
163  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=(op1(MA(m1,m2,m3)))
164 #End
165  end do
166  end do
167  end do
168 #endMacro
169 
170 #beginMacro setCoeff9(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9)
171 setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9,)
172 #endMacro
173 #beginMacro setCoeff8(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8)
174 setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,,)
175 #endMacro
176 #beginMacro setCoeff7(c,e,coeff,op1,op2,op3,op4,op5,op6,op7)
177 setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,,,)
178 #endMacro
179 #beginMacro setCoeff6(c,e,coeff,op1,op2,op3,op4,op5,op6)
180 setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,,,,)
181 #endMacro
182 #beginMacro setCoeff5(c,e,coeff,op1,op2,op3,op4,op5)
183 setCoeff10(c,e,coeff,op1,op2,op3,op4,op5,,,,,)
184 #endMacro
185 #beginMacro setCoeff4(c,e,coeff,op1,op2,op3,op4)
186 setCoeff10(c,e,coeff,op1,op2,op3,op4,,,,,,)
187 #endMacro
188 #beginMacro setCoeff3(c,e,coeff,op1,op2,op3)
189 setCoeff10(c,e,coeff,op1,op2,op3,,,,,,,)
190 #endMacro
191 #beginMacro setCoeff2(c,e,coeff,op1,op2)
192 setCoeff10(c,e,coeff,op1,op2,,,,,,,,)
193 #endMacro
194 #beginMacro setCoeff1(c,e,coeff,op1)
195 setCoeff10(c,e,coeff,op1,,,,,,,,,)
196 #endMacro
197 
198 
199 ! ======================================================================================================
200 ! Add up to 10 different local matrix operators to the global matrix coeff
201 ! in equation "e" and component "c"
202 !
203 ! (Leave final unused arguments empty)
204 ! ======================================================================================================
205 #beginMacro addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9,op10)
206  do m3=-halfWidth3,halfWidth3
207  do m2=-halfWidth,halfWidth
208  do m1=-halfWidth,halfWidth
209 #If #op10 ne ""
210  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
211  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
212  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))+(op9(MA(m1,m2,m3)))+(op10(MA(m1,m2,m3)))
213 #Elif #op9 ne ""
214  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
215  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
216  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))+(op9(MA(m1,m2,m3)))
217 #Elif #op8 ne ""
218  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
219  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
220  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))+(op8(MA(m1,m2,m3)))
221 #Elif #op7 ne ""
222  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
223  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
224  (op6(MA(m1,m2,m3)))+(op7(MA(m1,m2,m3)))
225 #Elif #op6 ne ""
226  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
227  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))+\
228  (op6(MA(m1,m2,m3)))
229 #Elif #op5 ne ""
230  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
231  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))+(op5(MA(m1,m2,m3)))
232 #Elif #op4 ne ""
233  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+\
234  (op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))+(op4(MA(m1,m2,m3)))
235 #Elif #op3 ne ""
236  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))+(op3(MA(m1,m2,m3)))
237 #Elif #op2 ne ""
238  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+(op1(MA(m1,m2,m3)))+(op2(MA(m1,m2,m3)))
239 #Else
240  coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)=coeff(MCE(m1,m2,m3,c,e),i1,i2,i3)+(op1(MA(m1,m2,m3)))
241 #End
242  end do
243  end do
244  end do
245 #endMacro
246 
247 #beginMacro addCoeff9(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9)
248 addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,op9,)
249 #endMacro
250 #beginMacro addCoeff8(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8)
251 addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,op8,,)
252 #endMacro
253 #beginMacro addCoeff7(c,e,coeff,op1,op2,op3,op4,op5,op6,op7)
254 addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,op7,,,)
255 #endMacro
256 #beginMacro addCoeff6(c,e,coeff,op1,op2,op3,op4,op5,op6)
257 addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,op6,,,,)
258 #endMacro
259 #beginMacro addCoeff5(c,e,coeff,op1,op2,op3,op4,op5)
260 addCoeff10(c,e,coeff,op1,op2,op3,op4,op5,,,,,)
261 #endMacro
262 #beginMacro addCoeff4(c,e,coeff,op1,op2,op3,op4)
263 addCoeff10(c,e,coeff,op1,op2,op3,op4,,,,,,)
264 #endMacro
265 #beginMacro addCoeff3(c,e,coeff,op1,op2,op3)
266 addCoeff10(c,e,coeff,op1,op2,op3,,,,,,,)
267 #endMacro
268 #beginMacro addCoeff2(c,e,coeff,op1,op2)
269 addCoeff10(c,e,coeff,op1,op2,,,,,,,,)
270 #endMacro
271 #beginMacro addCoeff1(c,e,coeff,op1)
272 addCoeff10(c,e,coeff,op1,,,,,,,,,)
273 #endMacro
274 
275 
276 
277 
278 
279 ! ======================================================================================================
280 ! Set the local matrix operator opCoeff to the global matrix coeff in equation "e" and component "c"
281 ! ======================================================================================================
282 #beginMacro setCoeff(c,e,coeff,opCoeff)
283 #If $DIM == 2
284  do m2=-halfWidth,halfWidth
285  do m1=-halfWidth,halfWidth
286  coeff(mce2(m1,m2,0,c,e),i1,i2,i3)=opCoeff(ma2(m1,m2,0))
287  end do
288  end do
289 #Else
290  do m3=-halfWidth,halfWidth
291  do m2=-halfWidth,halfWidth
292  do m1=-halfWidth,halfWidth
293  coeff(mce3(m1,m2,m3,c,e),i1,i2,i3)=opCoeff(ma3(m1,m2,m3))
294  end do
295  end do
296  end do
297 #End
298 #endMacro
299 
300 ! ==========================================================================================
301 ! Evaluate the Jacobian and its derivatives (parametric and spatial).
302 ! aj : prefix for the name of the resulting jacobian variables,
303 ! e.g. ajrx, ajsy, ajrxx, ajsxy, ...
304 ! MAXDER : number of derivatives to evaluate.
305 ! ==========================================================================================
306 #beginMacro opEvalJacobianDerivatives(aj,MAXDER)
307 
308 #If $GRIDTYPE eq "curvilinear"
309  ! this next call will define the jacobian and its derivatives (parameteric and spatial)
310  #peval evalJacobianDerivatives(rsxy,i1,i2,i3,aj,$DIM,$ORDER,MAXDER)
311 
312 #End
313 
314 #endMacro
315 
316 ! ==========================================================================================
317 ! Evaluate the parametric derivatives of u.
318 ! u : evaluate derivatives of this function.
319 ! uc : component to evaluate
320 ! uu : prefix for the name of the resulting derivatives, e.g. uur, uus, uurr, ...
321 ! MAXDER : number of derivatives to evaluate.
322 ! ==========================================================================================
323 #beginMacro opEvalParametricDerivative(u,uc,uu,MAXDER)
324 #If $GRIDTYPE eq "curvilinear"
325  #peval evalParametricDerivativesComponents1(u,i1,i2,i3,uc, uu,$DIM,$ORDER,MAXDER)
326 #Else
327  uu=u(i1,i2,i3,uc) ! in the rectangular case just eval the solution
328 #End
329 #endMacro
330 
331 
332 ! ==========================================================================================
333 ! Evaluate a derivative. (assumes parametric derivatives have already been evaluated)
334 ! DERIV : name of the derivative. One of
335 ! x,y,z,xx,xy,xz,...
336 ! u : evaluate derivatives of this function.
337 ! uc : component to evaluate
338 ! uu : prefix for the name of the resulting derivatives (same name used with opEvalParametricDerivative)
339 ! aj : prefix for the name of the jacobian variables.
340 ! ud : derivative is assigned to this variable.
341 ! ==========================================================================================
342 #beginMacro getOp(DERIV, u,uc,uu,aj,ud )
343 
344  #If $GRIDTYPE eq "curvilinear"
345  #peval getDuD ## DERIV ## $DIM(uu,aj,ud) ! Note: The perl variables are evaluated when the macro is USED.
346  #Else
347  #peval ud = u ## DERIV ## $ORDER(i1,i2,i3,uc)
348  #End
349 
350 #endMacro
351 
352 
353 ! ==========================================================================================
354 ! Form the local coefficient matrix for an operator. (assumes parametric derivatives have already been evaluated)
355 ! DERIV : name of the operator. One of
356 ! x,y,z,xx,xy,xz,...
357 ! laplacian, rr,ss,tt, rrrr,ssss,tttt,
358 ! r2Dissipation=rr+ss[+tt]
359 ! r4Dissipation=rrrr+ssss[+tttt]
360 ! coeff : fill in this (local) coefficient matrix
361 ! AJ : prefix for the name of the jacobian variables.
362 ! ==========================================================================================
363 #beginMacro getCoeff( DERIV, coeff, aj )
364 
365  #If $GRIDTYPE eq "curvilinear" || #DERIV eq "identity" || #DERIV eq "r2Dissipation" || #DERIV eq "r4Dissipation"
366  #peval DERIV ## CoeffOrder$ORDER\Dim$DIM(coeff,aj) ! trouble if ## appears after a perl variable
367  #Else
368  #peval DERIV ## CoeffOrder$ORDER\Dim$DIM\Rectangular(coeff,aj)
369  #End
370 
371 #endMacro
372 
373 
374 
375 ! ==============================================================================================================
376 ! Fill in the coefficients for a derivative
377 ! ==============================================================================================================
378 #beginMacro fillCoeff(DERIV,coeff,c,e)
379 beginLoops()
380 
381  ! Evaluate the jacobian derivatives used by the coefficient and forward derivatives:
382  ! opEvalJacobianDerivatives(MAXDER) : MAXDER = max number of derivatives to precompute.
383  opEvalJacobianDerivatives(aj,1)
384 
385  ! evaluate the coeff operators
386 
387  getCoeff(DERIV, xCoeff,aj)
388 
389  addCoeff(c,e,coeff,xCoeff)
390 
391 endLoops()
392 #endMacro
393 
394 ! ==============================================================================================================
395 ! Fill in the coefficients for a derivative *CONSERVATIVE*
396 ! DERIV : divScalarGrad
397 ! ==============================================================================================================
398 #beginMacro fillCoeffConservative(DERIV,s,coeff,c,e)
399 beginLoops()
400 
401  getConservativeCoeff( DERIV,s,xCoeff )
402 
403  addCoeff(c,e,coeff,xCoeff)
404 
405 endLoops()
406 #endMacro
407 
408 
409 
410 
411 !=======================================================================================
412 ! /Description: Return the equation number for given indices
413 ! /n (input): component number ( n=0,1,..,numberOfComponents-1 )
414 ! /i1,i2,i3 (input): grid indices
415 ! /return value : The equation number.
416 !\end{SparseRepInclude.tex}
417 !=======================================================================================
418 #defineMacro indexToEquation( n,i1,i2,i3 ) (n+1+ \
419  numberOfComponentsForCoefficients*(i1-equationNumberBase1+\
420  equationNumberLength1*(i2-equationNumberBase2+\
421  equationNumberLength2*(i3-equationNumberBase3))) + equationOffset)
422 
423 
424 !===============================================================================================
425 ! /Description:
426 ! Assign row and column numbers to entries in a sparse matrix.
427 ! This routine is normally only used for assign equation numbers on CompositeGrids
428 ! when the equationNumber belongs to a point on a different MappedGrid.
429 ! Rows and columns in the sparse matrix are numbered according to the values of
430 ! (n,I1,I2,I3)
431 ! where n is the component number and (I1,I2,I3) are the coordinate indicies on the grid.
432 ! The component number n runs from 0 to the numberOfComponentsForCoefficients-1 and is used
433 ! when solving a system of equations.
434 !
435 ! /m (input): assign row/column values for the m''th entry in the sparse matrix
436 ! /na,I1a,I2a,I3a (input): defines the row(s)
437 ! /equationNumber (input): defines an equation number
438 !
439 !\end{SparseRepInclude.tex}
440 !===============================================================================================
441 #beginMacro setEquationNumber(m, ni,i1,i2,i3, nj,j1,j2,j3 )
442  equationNumber(m,i1,i2,i3)=indexToEquation( nj,j1,j2,j3)
443 #endMacro
444 
445 
446 
447 !===============================================================================================
448 ! /Description:
449 ! Specify the classification for a set of Index values
450 !\end{SparseRepInclude.tex}
451 !===============================================================================================
452 #beginMacro setClassify(n,i1,i2,i3, type)
453  classify(i1,i2,i3,n)=type
454 #endMacro
455 
456 ! =======================================================================
457 ! Macro to zero out the matrix coefficients for equations e1,e1+1,..,e2
458 ! =======================================================================
459 #beginMacro zeroMatrixCoefficients( coeff,e1,e2, i1,i2,i3 )
460 do m=ce(0,e1),ce(0,e2+1)-1
461  coeff(m,i1,i2,i3)=0.
462 end do
463 #endMacro
464 
465 ! ===============================================================================================
466 ! Add an extrapolation equation to the matrix
467 ! Macro args:
468 ! coeff : coefficient matrix to fill in.
469 ! c,e : fill in equation e, extrapolate component c
470 ! i1,i2,i3 : marks the boundary point, ghost points will be assigned to assign
471 ! orderOfExtrap : order of extrapolation
472 ! ghost : point line to extrapolate (ghost=1,2,..)
473 ! ===============================================================================================
474 #beginMacro fillMatrixExtrapolation(coeff, c,e,i1,i2,i3,orderOfExtrap,ghost)
475  if( orderOfExtrap.lt.1 .or. orderOfExtrap.gt.maxOrderOfExtrapolation )then
476  stop 7734
477  end if
478  i1m=i1-is1*(ghost) ! ghost point
479  i2m=i2-is2*(ghost)
480  i3m=i3-is3*(ghost)
481  do m=0,orderOfExtrap
482  j1=i1m+is1*m ! m-th point moving inward from the ghost point (i1,i2,i3)
483  j2=i2m+is2*m
484  j3=i3m+is3*m
485  mm = ce(c,e)+m
486  coeff(mm,i1m,i2m,i3m)=extrapCoeff(m,orderOfExtrap) ! m=0,1,2,..
487 
488  setEquationNumber(mm, e,i1m,i2m,i3m, c,j1,j2,j3 ) !macro to set equationNumber
489  end do
490  setClassify(e,i1m,i2m,i3m, extrapolation) !macro to set classify
491 #endMacro
492 
493 
494 ! ===============================================================================================
495 ! Add a Neumann or mixed BC to the matrix
496 ! a0*I + a1*D_n
497 ! Macro args:
498 ! coeff : coefficient matrix to fill in.
499 ! c,e : fill in equation e, extrapolate component c
500 ! i1,i2,i3 : boundary point (will assign equations on the ghost point)
501 ! an(0:2) : holds the outward normal vector
502 ! a0,a1 : coefficients of the mixed BC
503 ! ===============================================================================================
504 #beginMacro fillMatrixNeumann(coeff, c,e, i1,i2,i3, an, a0,a1 )
505 
506  ! Evaluate the jacobian derivatives used by the coefficient and forward derivatives:
507  ! opEvalJacobianDerivatives(MAXDER) : MAXDER = max number of derivatives to precompute.
508  opEvalJacobianDerivatives(aj,0)
509 
510  ! evaluate the coeff operators
511  ! getCoeff(identity, iCoeff,aj)
512  getCoeff(x, xCoeff,aj)
513  getCoeff(y, yCoeff,aj)
514  getCoeff(identity, iCoeff,aj)
515  #If $DIM == 3
516  getCoeff(z, zCoeff,aj)
517  #End
518 
519  i1m=i1-is1 ! ghost point
520  i2m=i2-is2
521  i3m=i3-is3
522 
523  beginMatrixLoops(m1,m2,m3)
524  m=MA(m1,m2,m3)
525  mm=MCE(m1,m2,m3,c,e)
526 
527  #If $DIM == 2
528  coeff(mm,i1m,i2m,i3m)=a1*(an(0)*xCoeff(m)+an(1)*yCoeff(m))+a0*iCoeff(m)
529  #Else
530  coeff(mm,i1m,i2m,i3m)=a1*(an(0)*xCoeff(m)+an(1)*yCoeff(m)+an(2)*zCoeff(m))\
531  +a0*iCoeff(m)
532  #End
533 
534  ! The equation for pt (e,i1m,i2m,i3m) is centered on (c,i1,i2,i3):
535  setEquationNumber(mm, e,i1m,i2m,i3m, c,i1+m1,i2+m2,i3+m3 ) !macro to set equationNumber
536 
537  endMatrixLoops()
538  setClassify(e,i1m,i2m,i3m, ghost1) !macro to set classify
539 
540 #endMacro
541 
542 ! ===============================================================================================
543 ! Add a Vector Symmetry BC to the matrix
544 ! Macro args:
545 ! coeff : coefficient matrix to fill in.
546 ! cmpu,eqnu : fill in equations eqnu,...,eqnu+nd-1 and components cmpu,...,cmpu+nd-1
547 ! i1,i2,i3 : boundary point, will assign ghost point
548 ! NOTES:
549 ! The vector symmetry condition is the normal component is even, the tangential components are odd functions:
550 ! E1: n.u(-1) = - n.u(1)
551 ! E2: t.u(-1) = t.u(1)
552 ! or -> [ u(-1) - (n.u(-1))n ] = [ u(1) - (n.u(1))n ]
553 ! Combine equtions:
554 ! [n.u(-1) + n.u(1)]n + [ u(-1) - (n.u(-1))n ] - [ u(1) - (n.u(1))n ] = 0
555 ! or
556 ! u(-1) - u(1) + 2*(n.u(1))n = 0
557 ! ===============================================================================================
558 #beginMacro fillMatrixVectorSymmetry(coeff, cmpu,eqnu, i1,i2,i3, an )
559 
560  ! write(*,'(" VS: i1,i2=",2i2," normal=",2f5.2)') i1,i2,an(0),an(1)
561 
562  i1m=i1-is1 ! ghost point
563  i2m=i2-is2
564  i3m=i3-is3
565 
566  do m=0,ndc-1
567  coeff(m,i1m,i2m,i3m)=0. ! init all elements to zero
568  end do
569  mv(0)=0
570  mv(1)=0
571  mv(2)=0
572  do e=eqnu,eqnu+nd-1
573  c=cmpu+e-eqnu
574 
575  mv(axis)=2*side-1
576  mm=MCE(mv(0),mv(1),mv(2),c,e) ! matrix index for ghost pt
577  coeff(mm,i1m,i2m,i3m)=1. ! coeff of ghost point
578  setEquationNumber(mm, e,i1m,i2m,i3m, c,i1+mv(0),i2+mv(1),i3+mv(2) ) !macro to set equationNumber
579 
580  mv(axis)=-(2*side-1)
581  mm=MCE(mv(0),mv(1),mv(2),c,e) ! matrix index for first pt inside
582  coeff(mm,i1m,i2m,i3m)=-1.
583  setEquationNumber(mm, e,i1m,i2m,i3m, c,i1+mv(0),i2+mv(1),i3+mv(2) ) !macro to set equationNumber
584 
585  ! now add on the term: 2*(n.u(1))n
586  do c=cmpu,cmpu+nd-1
587  mm=MCE(mv(0),mv(1),mv(2),c,e) ! matrix index for first pt inside
588  ! write(*,'(" VS: mm,i1,i2,e,c=",5i2," n,n=",2f5.2)') mm,i1,i2,e,c,an(c-cmpu),an(e-eqnu)
589  coeff(mm,i1m,i2m,i3m)=coeff(mm,i1m,i2m,i3m) + 2.*an(c-cmpu)*an(e-eqnu)
590  setEquationNumber(mm, e,i1m,i2m,i3m, c,i1+mv(0),i2+mv(1),i3+mv(2) ) !macro to set equationNumber
591  end do
592 
593  setClassify(e,i1m,i2m,i3m, ghost1) !macro to set classify
594  end do
595 #endMacro
596 
597 ! ===============================================================================================
598 ! Return the normal vector (an(0),an(1),an(2)) for a point (i1,i2,i3) on a face (side,axis)
599 ! This macro does nothing on Cartesian grids.
600 ! ===============================================================================================
601 #beginMacro getNormalForCurvilinearGrid(side,axis,i1,i2,i3)
602  #If $GRIDTYPE == "curvilinear"
603  ! get the outward normal for curvilinear grids
604  an(0)=rsxy(i1,i2,i3,axis,0)
605  an(1)=rsxy(i1,i2,i3,axis,1)
606  #If $DIM == 2
607  anNorm = (2*side-1)/max( epsX, sqrt( an(0)**2 + an(1)**2 ) )
608  an(0)=an(0)*anNorm
609  an(1)=an(1)*anNorm
610  #Else
611  an(2)=rsxy(i1,i2,i3,axis,2)
612  anNorm = (2*side-1)/max( epsX, sqrt( an(0)**2 + an(1)**2 + an(2)**2 ) )
613  an(0)=an(0)*anNorm
614  an(1)=an(1)*anNorm
615  an(2)=an(2)*anNorm
616  #End
617  #End
618 #endMacro
619 
620 
621 ! *************************************************************************************
622 
623 
624 ! ==============================================================================================================
625 ! Fill in the coefficients for an equation
626 ! ==============================================================================================================
627 #beginMacro fillCoeff()
628 #If $DIM == 2
629  #defineMacro MCE(m1,m2,m3,c,e) mce2(m1,m2,m3,c,e)
630  #defineMacro MA(m1,m2,m3) ma2(m1,m2,m3)
631 #Else
632  #defineMacro MCE(m1,m2,m3,c,e) mce3(m1,m2,m3,c,e)
633  #defineMacro MA(m1,m2,m3) ma3(m1,m2,m3)
634 #End
635 
636  ! the next macro is defined in insImpINS.bf, or insImpVP.bf, etc.
637  fillCoeffPDE()
638 
639 #endMacro
640 
641 ! ======================================================================
642 ! Macro to fill in the matrix with BC-s
643 ! OR evaluate the residual on the boundary
644 !
645 ! Implicit parameters:
646 ! $DIM : 2 or 3
647 ! $GRIDTYPE : "rectangular" or "curvilinear"
648 ! ======================================================================
649 #beginMacro fillMatrixBoundaryConditions()
650 if( fillCoefficients.eq.1 .or. evalResidualForBoundaryConditions.eq.1 )then
651 
652 #If $DIM == 2
653  #defineMacro MCE(m1,m2,m3,c,e) mce2(m1,m2,m3,c,e)
654  #defineMacro MA(m1,m2,m3) ma2(m1,m2,m3)
655 #Else
656  #defineMacro MCE(m1,m2,m3,c,e) mce3(m1,m2,m3,c,e)
657  #defineMacro MA(m1,m2,m3) ma3(m1,m2,m3)
658 #End
659 
660  indexRange(0,0)=n1a
661  indexRange(1,0)=n1b
662  indexRange(0,1)=n2a
663  indexRange(1,1)=n2b
664  indexRange(0,2)=n3a
665  indexRange(1,2)=n3b
666 
667  do axis=0,nd-1
668  do side=0,1
669  is1=0
670  is2=0
671  is3=0
672  if( axis.eq.0 )then
673  is1=1-2*side
674  n1a=indexRange(side,axis)
675  n1b=n1a
676  else if( axis.eq.1 )then
677  is2=1-2*side
678  n2a=indexRange(side,axis)
679  n2b=n2a
680  else
681  is3=1-2*side
682  n3a=indexRange(side,axis)
683  n3b=n3a
684  end if
685 
686  #If $GRIDTYPE == "rectangular"
687  ! define the outward normal for Cartesian grids
688  an(0)=0.
689  an(1)=0.
690  an(2)=0.
691  an(axis)=2*side-1
692  #End
693 
694  ! the next macro is defined in insImpINS.bf, or insImpVP.bf, etc.
695  fillMatrixBoundaryConditionsPDE()
696  getBoundaryResidualPDE()
697 
698 
699 !* if( side.eq.0 .and. axis.eq.0 )then
700 !* fillMatrixBoundaryConditionsOnAFaceINS(normal00)
701 !* getBoundaryResidualINS(normal00)
702 !* else if( side.eq.1 .and. axis.eq.0 )then
703 !* fillMatrixBoundaryConditionsOnAFaceINS(normal10)
704 !* getBoundaryResidualINS(normal10)
705 !* else if( side.eq.0 .and. axis.eq.1 )then
706 !* fillMatrixBoundaryConditionsOnAFaceINS(normal01)
707 !* getBoundaryResidualINS(normal10)
708 !* else if( side.eq.1 .and. axis.eq.1 )then
709 !* fillMatrixBoundaryConditionsOnAFaceINS(normal11)
710 !* getBoundaryResidualINS(normal11)
711 !* else if( side.eq.0 .and. axis.eq.2 )then
712 !* fillMatrixBoundaryConditionsOnAFaceINS(normal02)
713 !* getBoundaryResidualINS(normal02)
714 !* else if( side.eq.1 .and. axis.eq.2 )then
715 !* fillMatrixBoundaryConditionsOnAFaceINS(normal12)
716 !* getBoundaryResidualINS(normal12)
717 !* else
718 !* ! stop 2077
719 !* end if
720 
721  ! reset values
722  if( axis.eq.0 )then
723  n1a=indexRange(0,axis)
724  n1b=indexRange(1,axis)
725  else if( axis.eq.1 )then
726  n2a=indexRange(0,axis)
727  n2b=indexRange(1,axis)
728  else
729  n3a=indexRange(0,axis)
730  n3b=indexRange(1,axis)
731  end if
732 
733  end do ! side
734  end do ! axis
735 
736 end if
737 #endMacro
738 
739 ! ===================================================================================
740 ! Macro to evaluate the RHS
741 ! ===================================================================================
742 #beginMacro assignRHS(PDE)
743  ! the next macro is defined in insImpINS.bf, or insImpVP.bf, etc.
744  assignRHSPDE()
745 #endMacro
746 
747 
748 ! ==============================================================================================================
749 ! Compute the coefficient of the artificial dissipation
750 ! ==============================================================================================================
751 #beginMacro getArtificialDissipationCoeff(adCoeff, u0x,u0y,u0z, v0x,v0y,v0z, w0x,w0y,w0z )
752 
753  #If $DIM == 2
754  adCoeff = ad21 + cd22*( abs(u0x)+abs(u0y) + abs(v0x)+abs(v0y) )
755  #Else
756  adCoeff = ad21 + cd22*( abs(u0x)+abs(u0y)+abs(u0z) + abs(v0x)+abs(v0y)+abs(v0z) + abs(w0x)+abs(w0y)+abs(w0z) )
757  #End
758 
759 #endMacro
760 
761 ! Second order dissipation operators:
762 #defineMacro uDiss22(u,uc) (u(i1-1,i2,i3,uc)+u(i1,i2-1,i3,uc)+u(i1+1,i2,i3,uc)+u(i1,i2+1,i3,uc)-4.*u(i1,i2,i3,uc))
763 #defineMacro uDiss23(u,uc) (u(i1-1,i2,i3,uc)+u(i1,i2-1,i3,uc)+u(i1,i2,i3-1,uc)+u(i1+1,i2,i3,uc)+u(i1,i2+1,i3,uc)+u(i1,i2,i3+1,uc)-6.*u(i1,i2,i3,uc))
764 
765 
766 #beginMacro INSIMP(NAME,PDE)
767  subroutine NAME(nd,nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b,\
768  mask,xy,rsxy,radiusInverse, u, ndc, coeff, fe,fi,ul, gv,gvl,dw, ndMatProp,matIndex,matValpc,matVal, bc, \
769  boundaryCondition, ndbcd1a,ndbcd1b,ndbcd2a,ndbcd2b,ndbcd3a,ndbcd3b,ndbcd4a,ndbcd4b,bcData,\
770  nde, equationNumber, classify, \
771  nr1a,nr1b,nr2a,nr2b,nr3a,nr3b, \
772  ipar, rpar, pdb, ierr )
773 c======================================================================
774 c
775 c Incompressible Navier Stokes IMPlicit
776 c -------------------------------------
777 c
778 c 1. Build the coefficient matrix for implicit methods
779 c 2. Evaluate the right-hand-side and residual
780 c
781 c nd : number of space dimensions
782 c nd1a,nd1b,nd2a,nd2b,nd3a,nd3b : array dimensions
783 c
784 c mask :
785 c xy :
786 c rsxy :
787 c coeff(m,i1,i2,i3) : array holding the matrix coefficients
788 c u : holds the current solution, used to form the coeff matrix.
789 c fe : holds the explicit part when evaluating the RHS
790 c fi : holds the implicit part when evaluating the RHS
791 c ul : holds the linearized solution, used when evaluating the linearized operator and RHS
792 c gv : gridVelocity for moving grids
793 c dw : distance to the wall for some turbulence models
794 c
795 c======================================================================
796  implicit none
797  integer nd, ndc, n1a,n1b,n2a,n2b,n3a,n3b,
798  & nd1a,nd1b,nd2a,nd2b,nd3a,nd3b,nd4a,nd4b
799  integer nde,nr1a,nr1b,nr2a,nr2b,nr3a,nr3b
800 
801  real u(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
802 
803  real coeff(0:ndc-1,nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
804 
805  real fe(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
806  real fi(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
807 
808  real ul(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,nd4a:nd4b)
809 
810  real gv(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
811  real gvl(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
812  real dw(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
813  real xy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1)
814  real rsxy(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:nd-1,0:nd-1)
815  real radiusInverse(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
816 
817  integer mask(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
818  integer bc(0:1,0:2),boundaryCondition(0:1,0:2),indexRange(0:1,0:2),ierr
819 
820  integer ndbcd1a,ndbcd1b,ndbcd2a,ndbcd2b,ndbcd3a,ndbcd3b,ndbcd4a,ndbcd4b
821  real bcData(ndbcd1a:ndbcd1b,ndbcd2a:ndbcd2b,ndbcd3a:ndbcd3b,ndbcd4a:ndbcd4b)
822 
823  integer equationNumber(0:nde-1,nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
824  integer classify(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:*)
825 
826  integer ipar(0:*)
827  real rpar(0:*)
828 
829  double precision pdb ! pointer to data base
830 
831  ! -- arrays for variable material properties --
832  integer constantMaterialProperties
833  integer piecewiseConstantMaterialProperties
834  integer variableMaterialProperties
835  parameter( constantMaterialProperties=0,\
836  piecewiseConstantMaterialProperties=1,\
837  variableMaterialProperties=2 )
838  integer materialFormat,ndMatProp
839  integer matIndex(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b)
840  real matValpc(0:ndMatProp-1,0:*)
841  real matVal(nd1a:nd1b,nd2a:nd2b,nd3a:nd3b,0:*)
842 
843 c ---- local variables -----
844  integer c,e,i1,i2,i3,m1,m2,m3,j1,j2,j3,ghostLine,n,i1m,i2m,i3m,i1p,i2p,i3p,ndu
845  integer side,axis,is1,is2,is3,mm,eqnTemp,debug,ntdc
846  integer kd,kd3,orderOfAccuracy,gridIsMoving,orderOfExtrap,orderOfExtrapolation,orderOfExtrapolationForOutflow
847  integer numberOfComponentsForCoefficients,stencilSize
848  integer gridIsImplicit,implicitOption,implicitMethod,
849  & isAxisymmetric,use2ndOrderAD,use4thOrderAD,fillCoefficientsScalarSystem
850  integer pc,uc,vc,wc,sc,nc,kc,ec,tc,qc,vsc,rc,grid,m,advectPassiveScalar
851  real nu,dt,nuPassiveScalar,adcPassiveScalar
852  real gravity(0:2), thermalExpansivity, adcBoussinesq,kThermal
853  real dxi,dyi,dzi,dri,dsi,dti,dr2i,ds2i,dt2i
854  real ad21,ad22,ad41,ad42,cd22,cd42,adc,adCoeff,adCoeffl
855  real ad21n,ad22n,ad41n,ad42n,cd22n,cd42n
856  real yy,yEps, epsX
857  real an(0:2),anNorm, advectionCoefficient
858  integer checkForInflowAtOutFlow, outflowOption
859  integer ok,getInt,getReal
860  integer mv(0:2)
861  integer fillCoefficients,evalRightHandSide,evalResidual,evalResidualForBoundaryConditions
862 
863  real ugv,vgv,wgv, ugvl,vgvl,wgvl
864 
865  integer equationNumberBase1,equationNumberLength1,equationNumberBase2,equationNumberLength2,\
866  equationNumberBase3,equationNumberLength3,equationOffset
867 
868  integer gridType
869  integer rectangular,curvilinear
870  parameter( rectangular=0, curvilinear=1 )
871 
872  integer turbulenceModel,noTurbulenceModel
873  integer baldwinLomax,spalartAllmaras,kEpsilon,kOmega,largeEddySimulation
874  parameter (noTurbulenceModel=0,baldwinLomax=1,kEpsilon=2,kOmega=3,spalartAllmaras=4,largeEddySimulation=5 )
875 
876  integer computeAllTerms,
877  & doNotComputeImplicitTerms,
878  & computeImplicitTermsSeparately,
879  & computeAllWithWeightedImplicit
880 
881  parameter( computeAllTerms=0,
882  & doNotComputeImplicitTerms=1,
883  & computeImplicitTermsSeparately=2,
884  & computeAllWithWeightedImplicit=3 )
885 
886  ! *new*
887  real nuDt,aDt,bDt,nuImp,aExp,aImp,bImp,tImp,kDt,teDt
888  real ulterm,vlterm,wlterm,qlterm
889  real implicitFactor
890  integer nonlinearTermsAreImplicit,evalLinearizedDerivatives
891 
892  integer implicitVariation
893  integer implicitViscous, implicitAdvectionAndViscous, implicitFullLinearized
894  parameter( implicitViscous=0,
895  & implicitAdvectionAndViscous=1,
896  & implicitFullLinearized=2 )
897 
898  integer pdeModel,standardModel,BoussinesqModel,viscoPlasticModel,twoPhaseFlowModel
899  parameter( standardModel=0,BoussinesqModel=1,viscoPlasticModel=2,twoPhaseFlowModel=3 )
900 
901  integer \
902  noSlipWall, \
903  inflowWithVelocityGiven, \
904  slipWall, \
905  outflow, \
906  convectiveOutflow, \
907  tractionFree, \
908  inflowWithPandTV, \
909  dirichletBoundaryCondition, \
910  symmetry, \
911  axisymmetric, \
912  interfaceBoundaryCondition,\
913  neumannBoundaryCondition
914  parameter( noSlipWall=1,inflowWithVelocityGiven=2,slipWall=4, \
915  outflow=5,convectiveOutflow=14,tractionFree=15, \
916  inflowWithPandTV=3, \
917  dirichletBoundaryCondition=12, \
918  symmetry=11,axisymmetric=13,interfaceBoundaryCondition=17,neumannBoundaryCondition=101 )
919 
920  ! classifyTypes from SparseRep.h
921 
922  integer interior,boundary,ghost1,ghost2,ghost3,ghost4,interpolation,periodic,extrapolation,unused
923  parameter(interior=1,\
924  boundary=2,\
925  ghost1=3,\
926  ghost2=4,\
927  ghost3=5,\
928  ghost4=6,\
929  interpolation=-1,\
930  periodic=-2,\
931  extrapolation=-3,\
932  unused=0 )
933 
934  ! the following define which component to fill for a scalar coeff problem
935  integer fillCoeffU,fillCoeffV,fillCoeffW,fillCoeffT
936  parameter( fillCoeffU=1,fillCoeffV=2,fillCoeffW=3,fillCoeffT=4 )
937 
938  real rx,ry,rz,sx,sy,sz,tx,ty,tz
939  real aj
940  real dr(0:2), dx(0:2)
941 
942  integer ma2(-5:5,-5:5,0:0), ma3(-5:5,-5:5,-5:5)
943  real nDotL(0:2)
944  real delta(0:5,0:5)
945  integer ncc,width,halfWidth,halfWidth3
946 
947  integer maxWidth,maxWidthDim
948  parameter( maxWidth=5,maxWidthDim=maxWidth*maxWidth*maxWidth )
949  real iCoeff(0:maxWidthDim-1)
950  real xCoeff(0:maxWidthDim-1)
951  real yCoeff(0:maxWidthDim-1)
952  real zCoeff(0:maxWidthDim-1)
953  real lapCoeff(0:maxWidthDim-1)
954  real divSGradCoeff(0:maxWidthDim-1)
955  real dissCoeff(0:maxWidthDim-1)
956 
957  real p0x,p0y,p0z
958  real ulx,uly,ulz, vlx,vly,vlz, wlx,wly,wlz, plx,ply,plz, qlx,qly,qlz
959  real u0x,u0y,u0z,u0xx,u0xy,u0xz,u0yy,u0yz,u0zz
960  real v0x,v0y,v0z,v0xx,v0xy,v0xz,v0yy,v0yz,v0zz
961  real w0x,w0y,w0z,w0xx,w0xy,w0xz,w0yy,w0yz,w0zz
962  real q0x,q0y,q0z,q0xx,q0xy,q0xz,q0yy,q0yz,q0zz
963  integer eqn,eqnu,eqnv,eqnw,eqnq,eqnk,eqne
964  integer cmp,cmpu,cmpv,cmpw,cmpq,cmpk,cmpe
965  integer mce2,mce3,ce
966  ! temp variables for coeff macros:
967  real cur,cus,cut,curr,curs,curt,cuss,cust,cutt
968 
969  ! --- visco plastic variables ---
970  real dr0i,dr1i,dr0dr1,dr0dr2,dr1dr2
971  real divNuGradu,divNuGradv,divNuGradw,divNuGradul,divNuGradvl,divNuGradwl
972  real dxvsqi(0:2),dtImp
973 
974  ! This include file (created in insImpINS.bf or insImpVP, ...) declares variables needed by each version
975  declareInsImpTemporaryVariables()
976 
977  integer maxOrderOfExtrapolation
978  parameter( maxOrderOfExtrapolation=9 )
979  real extrapCoeff(0:maxOrderOfExtrapolation,1:maxOrderOfExtrapolation)
980 
981  integer numberOfComponents
982  real mixedRHS,mixedCoeff,mixedNormalCoeff,a0,a1
983 
984  real rhopc,rhov, Cppc, Cpv, thermalKpc, thermalKv, Kx, Ky, Kz, Kr, Ks, Kt
985 
986  ! --- begin statement functions
987  rx(i1,i2,i3)=rsxy(i1,i2,i3,0,0)
988  ry(i1,i2,i3)=rsxy(i1,i2,i3,0,1)
989  rz(i1,i2,i3)=rsxy(i1,i2,i3,0,2)
990  sx(i1,i2,i3)=rsxy(i1,i2,i3,1,0)
991  sy(i1,i2,i3)=rsxy(i1,i2,i3,1,1)
992  sz(i1,i2,i3)=rsxy(i1,i2,i3,1,2)
993  tx(i1,i2,i3)=rsxy(i1,i2,i3,2,0)
994  ty(i1,i2,i3)=rsxy(i1,i2,i3,2,1)
995  tz(i1,i2,i3)=rsxy(i1,i2,i3,2,2)
996 
997  ! c=component, e=equation
998  ! ncc = number of components for coefficients
999  mce2(m1,m2,m3,c,e) = ma2(m1,m2,m3) + stencilSize*(c+ncc*e)
1000  mce3(m1,m2,m3,c,e) = ma3(m1,m2,m3) + stencilSize*(c+ncc*e)
1001  ce(c,e) = stencilSize*(c+ncc*e)
1002 
1003  ! statement functions to access coefficients of mixed-boundary conditions
1004  mixedRHS(c,side,axis,grid) =bcData(c+numberOfComponents*(0),side,axis,grid)
1005  mixedCoeff(c,side,axis,grid) =bcData(c+numberOfComponents*(1),side,axis,grid)
1006  mixedNormalCoeff(c,side,axis,grid) =bcData(c+numberOfComponents*(2),side,axis,grid)
1007 
1008 
1009  ! -- statement functions for variable material properties
1010  ! (rho,Cp,k) for materialFormat=piecewiseConstantMaterialProperties
1011  rhopc(i1,i2,i3) = matValpc( 0, matIndex(i1,i2,i3))
1012  Cppc(i1,i2,i3) = matValpc( 1, matIndex(i1,i2,i3))
1013  thermalKpc(i1,i2,i3) = matValpc( 2, matIndex(i1,i2,i3))
1014 
1015  ! (rho,Cp,k) for materialFormat=variableMaterialProperties
1016  rhov(i1,i2,i3) = matVal(i1,i2,i3,0)
1017  Cpv(i1,i2,i3) = matVal(i1,i2,i3,1)
1018  thermalKv(i1,i2,i3) = matVal(i1,i2,i3,2)
1019 
1020  ! --- end statement functions
1021 
1022  data extrapCoeff/ \
1023  1.,-1.,0.,0.,0.,0.,0.,0.,0.,0., \
1024  1.,-2.,1.,0.,0.,0.,0.,0.,0.,0., \
1025  1.,-3.,3.,-1.,0.,0.,0.,0.,0.,0., \
1026  1.,-4.,6.,-4.,1.,0.,0.,0.,0.,0., \
1027  1.,-5.,10.,-10.,5.,-1.,0.,0.,0.,0., \
1028  1.,-6.,15.,-20.,15.,-6.,1.,0.,0.,0., \
1029  1.,-7.,21.,-35.,35.,-21.,7.,-1.,0.,0., \
1030  1.,-8.,28.,-56.,70.,-56.,28.,-8.,1.,0., \
1031  1.,-9.,36.,-84.,126.,-126.,84.,-36.,9.,-1./
1032  ierr=0
1033  ! write(*,'("Inside insdt: gridType=",i2)') gridType
1034 
1035  n1a =ipar(0)
1036  n1b =ipar(1)
1037  n2a =ipar(2)
1038  n2b =ipar(3)
1039  n3a =ipar(4)
1040  n3b =ipar(5)
1041 
1042  pc =ipar(6)
1043  uc =ipar(7)
1044  vc =ipar(8)
1045  wc =ipar(9)
1046  nc =ipar(10)
1047  sc =ipar(11)
1048  tc =ipar(12)
1049  grid =ipar(13)
1050  orderOfAccuracy =ipar(14)
1051  gridIsMoving =ipar(15) ! *************
1052 
1053  implicitVariation =ipar(16) ! **new**
1054 
1055  fillCoefficients =ipar(17) ! new
1056  evalRightHandSide =ipar(18) ! new
1057 
1058  gridIsImplicit =ipar(19)
1059  implicitMethod =ipar(20)
1060  implicitOption =ipar(21)
1061  isAxisymmetric =ipar(22)
1062  use2ndOrderAD =ipar(23)
1063  use4thOrderAD =ipar(24)
1064  advectPassiveScalar=ipar(25)
1065  gridType =ipar(26)
1066  turbulenceModel =ipar(27)
1067  pdeModel =ipar(28)
1068  numberOfComponentsForCoefficients =ipar(29) ! number of components for coefficients
1069  stencilSize =ipar(30)
1070 
1071  equationOffset = ipar(31)
1072  equationNumberBase1 = ipar(32)
1073  equationNumberLength1 = ipar(33)
1074  equationNumberBase2 = ipar(34)
1075  equationNumberLength2 = ipar(35)
1076  equationNumberBase3 = ipar(36)
1077  equationNumberLength3 = ipar(37)
1078 
1079  indexRange(0,0) =ipar(38)
1080  indexRange(1,0) =ipar(39)
1081  indexRange(0,1) =ipar(40)
1082  indexRange(1,1) =ipar(41)
1083  indexRange(0,2) =ipar(42)
1084  indexRange(1,2) =ipar(43)
1085 
1086  orderOfExtrapolation=ipar(44)
1087  orderOfExtrapolationForOutflow=ipar(45)
1088 
1089  evalResidual = ipar(46)
1090  evalResidualForBoundaryConditions=ipar(47)
1091  debug = ipar(48)
1092  numberOfComponents= ipar(49)
1093 
1094  rc =ipar(50)
1095  materialFormat =ipar(51)
1096 
1097  dr(0) =rpar(0)
1098  dr(1) =rpar(1)
1099  dr(2) =rpar(2)
1100  dx(0) =rpar(3)
1101  dx(1) =rpar(4)
1102  dx(2) =rpar(5)
1103 
1104  dt =rpar(6) ! **new**
1105  implicitFactor =rpar(7) ! **new**
1106 
1107  nu =rpar(8)
1108  ad21 =rpar(9)
1109  ad22 =rpar(10)
1110  ad41 =rpar(11)
1111  ad42 =rpar(12)
1112  nuPassiveScalar =rpar(13)
1113  adcPassiveScalar =rpar(14)
1114  ad21n =rpar(15)
1115  ad22n =rpar(16)
1116  ad41n =rpar(17)
1117  ad42n =rpar(18)
1118  yEps =rpar(19) ! for axisymmetric
1119 
1120  gravity(0) =rpar(20)
1121  gravity(1) =rpar(21)
1122  gravity(2) =rpar(22)
1123  thermalExpansivity=rpar(23)
1124  adcBoussinesq =rpar(24) ! coefficient of artificial diffusion for Boussinesq T equation
1125  kThermal =rpar(25)
1126 
1127  epsX = 1.e-30 ! epsilon used to avoid division by zero in the normal computation -- should be REAL_MIN*100 ??
1128 
1129  ncc=numberOfComponentsForCoefficients ! number of components for coefficients
1130 
1131  ok = getInt(pdb,'checkForInflowAtOutFlow',checkForInflowAtOutFlow)
1132  if( ok.eq.0 )then
1133  write(*,'("*** NAME:ERROR: checkForInflowAtOutFlow NOT FOUND")')
1134  else
1135  if( debug.gt.4 )then
1136  write(*,'("*** NAME:checkForInflowAtOutFlow=",i6)') checkForInflowAtOutFlow
1137  end if
1138  end if
1139  ok = getInt(pdb,'outflowOption',outflowOption)
1140  if( ok.eq.0 )then
1141  write(*,'("*** NAME:ERROR: outflowOption NOT FOUND")')
1142  else
1143  if( debug.gt.4 )then
1144  write(*,'("*** NAME:outflowOption=",i6)') outflowOption
1145  end if
1146  end if
1147 
1148  ok = getReal(pdb,'advectionCoefficient',advectionCoefficient)
1149  if( ok.eq.0 )then
1150  write(*,'("*** NAME:ERROR: advectionCoefficient NOT FOUND")')
1151  else
1152  if( debug.gt.4 )then
1153  write(*,'("*** NAME:advectionCoefficient=",f5.2)') advectionCoefficient
1154  end if
1155  end if
1156 
1157 
1158  ok = getInt(pdb,'vsc',vsc)
1159  if( ok.eq.0 )then
1160  write(*,'("*** NAME:ERROR: vsc NOT FOUND")')
1161  end if
1162 
1163  ntdc = nd ! number of time dependent components
1164  if( pdeModel.eq.BoussinesqModel .or. pdeModel.eq.viscoPlasticModel )then
1165  ntdc = ntdc+1 ! include Temperature
1166  end if
1167 
1168  kc=nc
1169  ec=kc+1
1170 
1171  ok = getInt(pdb,'fillCoefficientsScalarSystem',fillCoefficientsScalarSystem)
1172  if( ok.eq.0 )then
1173  write(*,'("*** NAME:ERROR: fillCoefficientsScalarSystem NOT FOUND")')
1174  else
1175  if( debug.gt.4 )then
1176  write(*,'("*** NAME:fillCoefficientsScalarSystem=",i6)') fillCoefficientsScalarSystem
1177  end if
1178  end if
1179 
1180  if( orderOfAccuracy.ne.2 .and. orderOfAccuracy.ne.4 )then
1181  write(*,'("NAME:ERROR orderOfAccuracy=",i6)') orderOfAccuracy
1182  stop 1
1183  end if
1184  if( gridType.ne.rectangular .and. gridType.ne.curvilinear )then
1185  write(*,'("NAME:ERROR gridType=",i6)') gridType
1186  stop 2
1187  end if
1188  if( uc.lt.0 .or. vc.lt.0 .or. (nd.eq.3 .and. wc.lt.0) )then
1189  write(*,'("NAME:ERROR uc,vc,ws=",3i6)') uc,vc,wc
1190  stop 4
1191  end if
1192 
1193 c write(*,'("NAME: turbulenceModel=",2i6)') turbulenceModel
1194 c write(*,'("NAME: nd,uc,vc,wc,kc=",2i6)') nd,uc,vc,wc,kc
1195 
1196  if( turbulenceModel.eq.kEpsilon .and. (kc.lt.uc+nd .or. kc.gt.1000) )then
1197  write(*,'("NAME:ERROR in kc: nd,uc,vc,wc,kc=",2i6)') nd,uc,vc,wc,kc
1198  stop 5
1199  end if
1200 
1201 
1202 ! *****************************************************************************************************8
1203  ! ** it did not affect performance to use an array to index coeff ***
1204 
1205  width = orderOfAccuracy+1 ! width of the MATRIX stencil *********** fix this **********
1206  halfWidth = (width-1)/2
1207 
1208  if( nd.eq.2 )then
1209  halfWidth3=0
1210  do i2=-halfWidth,halfWidth
1211  do i1=-halfWidth,halfWidth
1212  ma2(i1,i2,0)=i1+halfWidth+width*(i2+halfWidth)
1213  end do
1214  end do
1215  else if( nd.eq.3 )then
1216  halfWidth3=halfWidth
1217  do i3=-halfWidth,halfWidth
1218  do i2=-halfWidth,halfWidth
1219  do i1=-halfWidth,halfWidth
1220  ma3(i1,i2,i3)=i1+halfWidth+width*(i2+halfWidth+width*(i3+halfWidth))
1221  end do
1222  end do
1223  end do
1224  end if
1225 
1226  do e=0,5
1227  do c=0,5
1228  if( e.eq.c )then
1229  delta(c,e)=1.
1230  else
1231  delta(c,e)=0.
1232  end if
1233  end do
1234  end do
1235 
1236  do m=0,2
1237  dxvsqi(m)=1./(dx(m)**2)
1238  end do
1239 
1240  if( use2ndOrderAD.ne.0 )then
1241  if( debug .gt.3 )then
1242  write(*,'(" NAME: INFO: 2nd order art-diss is on: ad21,ad22=",2(e9.2,1x))') ad21,ad22
1243  end if
1244  else
1245  ad21=0.
1246  ad22=0.
1247  end if
1248  cd22=ad22/(nd**2) ! for 2nd-order artificial dissipation
1249 
1250  if( pdeModel.eq.BoussinesqModel )then
1251  if( debug .gt.3 )then
1252  write(*,'(" NAME: INFO: Boussinesq: kThermal,adcBoussinesq=",2(e9.2,1x))') kThermal,adcBoussinesq
1253  end if
1254  end if
1255 
1256 
1257  ! *** define constants in the implicit operator
1258  !
1259 
1260  cmpu=uc-uc ! component number for u in the matrix
1261  cmpv=vc-uc
1262  cmpw=wc-uc
1263  qc=tc ! another name for tc
1264  cmpq=tc-uc ! temperature
1265  cmpk=kc-uc ! k
1266  cmpe=ec-uc ! eps
1267 
1268  ndu=nd ! number of velocity components in the matrix
1269  if( fillCoefficientsScalarSystem.ge.fillCoeffU .and. fillCoefficientsScalarSystem.le.fillCoeffW )then
1270  ! we are filling a scalar matrix for one component of the velocity
1271  ndu=1
1272  cmpv=cmpu
1273  cmpw=cmpu
1274  cmpq=cmpu
1275  else if( fillCoefficientsScalarSystem.eq.fillCoeffT )then
1276  ! we are filling a scalar matrix for T
1277  ndu=0
1278  cmpv=cmpu
1279  cmpw=cmpu
1280  cmpq=cmpu
1281  end if
1282 
1283  eqnu=cmpu ! equation number for u in the matrix
1284  eqnv=cmpv
1285  eqnw=cmpw
1286  eqnq=cmpq
1287  eqnk=cmpk
1288  eqne=cmpe
1289 
1290  dtImp=dt*implicitFactor
1291 
1292  nuDt= nu*dt*implicitFactor ! matrix: coefficient of Laplacian in the matrix
1293  kDt = kThermal*dt*implicitFactor ! for T equation
1294  aDt = dt*implicitFactor ! matrix: coefficient of advection terms in the matrix
1295  bDt = dt*implicitFactor ! matrix: coefficient of extra zero-order linearized terms such as u0.x
1296  teDt = dt*thermalExpansivity*implicitFactor ! matrix: coeff of buoyancy term
1297 
1298  nuImp=nu ! RHS : nu
1299  aImp=1. ! RHS: coeff of implicit advection terms
1300  bImp=1. ! RHS: coefficient of extra zero-order linearized terms such as u0.x
1301  ! nuImp=nu*(1.-implicitFactor) ! for RHS
1302  ! aImp=(1.-implicitFactor) ! RHS: coeff of implicit advection terms
1303  ! aExp=0. ! advection terms are implicit
1304 
1305  ! ** NOTE: we force the buoyancy terms to be explicit below **
1306  tImp=1. ! coefficient of the implicit buoyancy term (set to zero if not implicit)
1307  if( fillCoefficientsScalarSystem.ne.0 )then
1308  tImp=0. ! buoyancy term is not implicit when we solve scalar systems
1309  end if
1310 
1311 
1312  nonlinearTermsAreImplicit=0
1313  if( implicitVariation.eq.implicitViscous )then
1314  aDt=0. ! matrix: advection terms are NOT implicit
1315  bDt=0. ! matrix: zero-order linearized terms are NOT implicit
1316  ! aExp=1. ! advection terms are explicit
1317  aImp=0. ! RHS: coeff of implicit advection terms
1318  bImp=0.
1319  tImp=0. ! Boussinesq terms are explicit
1320  else if( implicitVariation.eq.implicitAdvectionAndViscous )then
1321  nonlinearTermsAreImplicit=1
1322  bDt=0.
1323  bImp=0.
1324  tImp=1.
1325  else if( implicitVariation.eq.implicitFullLinearized )then
1326  nonlinearTermsAreImplicit=1
1327  ! ++ bDt=0.
1328  ! ++ bImp=0.
1329  tImp=1.
1330  else
1331  write(*,'(" NAME: ERROR: unexpected implicitVariation=",i6)') implicitVariation
1332  stop 7200
1333  end if
1334 
1335  ! *** force the buoyancy term to be explicit : (we could add this as an option)
1336  ! if the main balance in the momentum is grad(p) = - alpha*g*T then we want both these
1337  ! terms to be explicit.
1338  teDt=0.
1339  tImp=0.
1340 
1341 
1342  ! The next macro is defined in insImpXX.bf and is used to look up parameters etc.
1343  initializePdeParameters()
1344 
1345 
1346  if( debug.gt.3 )then
1347  if( evalRightHandSide.eq.1 )then
1348  write(*,'("NAME: *EVAL RHS : implicitOption=",i2," (0=all,1=none,2=sep,3=all)")') implicitOption
1349  endif
1350  endif
1351  if( debug.gt.7 )then
1352  if( evalRightHandSide.eq.1 )then
1353  write(*,'("NAME: ******** EVAL RHS **********")')
1354  end if
1355  if( evalResidual.eq.1 )then
1356  write(*,'("NAME: ******** EVAL RESIDUAL **********")')
1357  end if
1358  if( fillCoefficients.eq.1 )then
1359  write(*,'("NAME: ******** FILL COEFF **********")')
1360  end if
1361  ! '
1362 
1363  write(*,'("NAME: pdeModel,nonLinear,impVar=",3i4)') pdeModel,nonlinearTermsAreImplicit,implicitVariation
1364  if( fillCoefficients.eq.1 )then
1365  write(*,'("NAME: stencilSize,ncc=",2i4)') stencilSize,ncc
1366  end if
1367  write(*,'("NAME: aDt,bDt,nuDt=",3e10.2)') aDt,bDt,nuDt
1368  write(*,'("NAME: implicitFactor,nuImp,aImp=",f4.2,1x,3e10.2)') implicitFactor,nuImp,aImp
1369  !'
1370  end if
1371 
1372  if( evalRightHandSide.eq.1 .and. (nonlinearTermsAreImplicit.eq.1 .or. use2ndOrderAD.ne.0) )then
1373  evalLinearizedDerivatives = 1
1374  else
1375  evalLinearizedDerivatives = 0
1376  end if
1377 
1378 
1379 
1380 ! Define operator parameters
1381 ! $DIM : number of spatial dimensions
1382 ! $ORDER : order of accuracy of an approximation
1383 ! $GRIDTYPE : rectangular or curvilinear
1384 ! $MATRIX_STENCIL_WIDTH : space in the global coeff matrix was allocated to hold this size stencil
1385 ! $STENCIL_WIDTH : stencil width of the local coeff-matrix (such as xCoeff, yCoeff, lapCoeff, ...)
1386 #perl $ORDER=2; $MATRIX_STENCIL_WIDTH=3; $STENCIL_WIDTH=3;
1387 
1388  if( nd.eq.2 .and. gridType.eq.rectangular .and. orderOfAccuracy.eq.2 )then
1389 #perl $DIM=2; $GRIDTYPE="rectangular";
1390 
1391  ! fill the coefficients:
1392  fillCoeff()
1393  ! fill matrix BCs
1394  fillMatrixBoundaryConditions()
1395  ! assign the RHS:
1396  assignRHS()
1397  else if( nd.eq.2 .and. gridType.eq.curvilinear .and. orderOfAccuracy.eq.2 )then
1398 #perl $DIM=2; $GRIDTYPE="curvilinear";
1399 
1400  ! fill the coefficients:
1401  fillCoeff()
1402  ! fill matrix BCs
1403  fillMatrixBoundaryConditions()
1404  ! assign the RHS:
1405  assignRHS()
1406 
1407  else if( nd.eq.3 .and. gridType.eq.rectangular .and. orderOfAccuracy.eq.2 )then
1408 #perl $DIM=3; $GRIDTYPE="rectangular";
1409 
1410  ! fill the coefficients:
1411  fillCoeff()
1412  ! fill matrix BCs
1413  fillMatrixBoundaryConditions()
1414  ! assign the RHS:
1415  assignRHS()
1416 
1417  else if( nd.eq.3 .and. gridType.eq.curvilinear .and. orderOfAccuracy.eq.2 )then
1418 #perl $DIM=3; $GRIDTYPE="curvilinear";
1419 
1420  ! fill the coefficients:
1421  fillCoeff()
1422  ! fill matrix BCs
1423  fillMatrixBoundaryConditions()
1424  ! assign the RHS:
1425  assignRHS()
1426 
1427  else
1428 
1429  if( orderOfAccuracy.ne.2 )then
1430  write(*,'("NAME: ERROR - not implemented for orderOfAccuracy=",i4)') orderOfAccuracy
1431  end if
1432 
1433 
1434  stop 9425
1435  end if
1436 
1437  return
1438  end
1439 
1440 #endMacro
1441 
Generated on Fri Jan 4 2013 10:19:22 for CG by   doxygen 1.8.3