CG  Version 25
lineSolveVP.h
Go to the documentation of this file.
1 c ************************************************************************
2 c Visco-Plastic Model:
3 c Macros to define the equations for the line solver
4 c
5 c This file is included by insLineSolveNew.bf
6 c
7 c In 2d the momentum equations are:
8 c
9 c u.t + u.grad(u) + p.x = Dx( 2*vp*u.x ) + Dy( vp*u.y ) + Dy( vp*v.x )
10 c v.t + u.grad(v) + p.y = Dx( vp*v.x ) + Dy( 2*vp*v.y ) + Dx( vp*u.y )
11 c
12 c ************************************************************************
13 
14 ! The next include file defines conservative approximations to coefficent matrices
15 #Include consCoeff.h
16 
17 
18 c ===================================================================================
19 c VP: incompressible Navier Stokes with a visco-plastic model
20 c *** rectangular grid version ***
21 c Macro arguments:
22 c dir: 0,1,2
23 c====================================================================================
24 #beginMacro fillEquationsRectangularGridVP(dir)
25 
26  if( use4thOrderAD.eq.1 )then
27  write(*,*) 'insLineSolve: 4th order diss not finished'
28  stop 7654
29  end if
30 
31  ! set default values for no 2nd order artificial diffusion:
32  cdm=0.
33  cdDiag=0.
34  cdp=0.
35 
36 if( nd.eq.2 )then
37 
38  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
39  defineDerivativeMacros(2,2,rectangular)
40 
41 
42  beginLoops()
43  if( mask(i1,i2,i3).gt.0 )then
44  if( use2ndOrderAD.eq.1 )then
45  defineArtificialDiffusionCoefficients(2,dir,R,)
46  end if
47 
48  ! Get the nonlinear viscosity at nearby points:
49  nuzmz=u(i1 ,i2-1,i3,vsc)
50  numzz=u(i1-1,i2 ,i3,vsc)
51  nuzzz=u(i1 ,i2 ,i3,vsc)
52  nupzz=u(i1+1,i2 ,i3,vsc)
53  nuzpz=u(i1 ,i2+1,i3,vsc)
54  ! getViscoPlasticViscosityCoefficient(nuzmz,i1 ,i2-1,i3,2,rectangular)
55  ! getViscoPlasticViscosityCoefficient(numzz,i1-1,i2 ,i3,2,rectangular)
56  ! getViscoPlasticViscosityCoefficient(nuzzz,i1 ,i2 ,i3,2,rectangular)
57  ! getViscoPlasticViscosityCoefficient(nupzz,i1+1,i2 ,i3,2,rectangular)
58  ! getViscoPlasticViscosityCoefficient(nuzpz,i1 ,i2+1,i3,2,rectangular)
59 
60  ! u.t + u.grad(u) + p.x = Dx( 2*nu*u.x ) + Dy( nu*u.y ) + Dy( nu*v.x )
61  ! v.t + u.grad(v) + p.y = Dx( nu*v.x ) + Dy( 2*nu*v.y ) + Dx( nu*u.y )
62 
63 
64  nu0ph = .5*( nupzz+nuzzz ) ! nu(i1+1/2,i2,i3)
65  nu0mh = .5*( nuzzz+numzz ) ! nu(i1-1/2,i2,i3)
66 
67  nu1ph = .5*( nuzpz+nuzzz ) ! nu(i1,i2+1/2,i3)
68  nu1mh = .5*( nuzzz+nuzmz ) ! nu(i1,i2-1/2,i3)
69 
70  if( computeMatrix.eq.1 )then
71  #If #dir == "0"
72  if( systemComponent.eq.uc )then ! solve for u
73  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir) -cdm - 2.*nu0mh*dxvsqi(0)
74  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + 2.*(nu0ph+nu0mh)*dxvsqi(0)+(nu1ph+nu1mh)*dxvsqi(1)
75  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir) -cdp - 2.*nu0ph*dxvsqi(0)
76  else if( systemComponent.eq.vc )then ! solve for v
77  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir) -cdm - nu0mh*dxvsqi(0)
78  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (nu0ph+nu0mh)*dxvsqi(0)+2.*(nu1ph+nu1mh)*dxvsqi(1)
79  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir) -cdp - nu0ph*dxvsqi(0)
80  else
81  stop 6139
82  end if
83  #Else
84  if( systemComponent.eq.uc )then ! solve for u
85  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir) -cdm - nu1mh*dxvsqi(1)
86  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + 2.*(nu0ph+nu0mh)*dxvsqi(0)+(nu1ph+nu1mh)*dxvsqi(1)
87  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir) -cdp - nu1ph*dxvsqi(1)
88  else if( systemComponent.eq.vc )then ! solve for v
89  am(i1,i2,i3)= -uu(uc+dir)*dxv2i(dir) -cdm - 2.*nu1mh*dxvsqi(1)
90  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (nu0ph+nu0mh)*dxvsqi(0)+2.*(nu1ph+nu1mh)*dxvsqi(1)
91  cm(i1,i2,i3)= uu(uc+dir)*dxv2i(dir) -cdp - 2.*nu1ph*dxvsqi(1)
92  else
93  stop 6139
94  end if
95  #End
96 
97  end if
98 
99 ! Here is the full residual:
100 ! f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu) + uu(uc)*dtScale/dt(i1,i2,i3) \
101 ! -u(i1,i2,i3,uc)*ux2(uc)-u(i1,i2,i3,vc)*uy2(uc)-ux2(pc)\
102 ! -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
103 ! +2.*(nu0ph*u(i1+1,i2,i3,uc) -(nu0ph+nu0mh)*u(i1,i2,i3,uc) + nu0mh*u(i1-1,i2,i3,uc))*dxvsqi(0)\
104 ! + (nu1ph*u(i1,i2+1,i3,uc) -(nu1ph+nu1mh)*u(i1,i2,i3,uc) + nu1mh*u(i1,i2-1,i3,uc))*dxvsqi(1)\
105 ! + ( nuzpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))\
106 ! -nuzmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dx(0)*dx(1))
107 !
108 ! f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv) + uu(vc)*dtScale/dt(i1,i2,i3) \
109 ! -u(i1,i2,i3,uc)*ux2(vc)-u(i1,i2,i3,vc)*uy2(vc)-uy2(pc)\
110 ! -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
111 ! + (nu0ph*u(i1+1,i2,i3,vc) -(nu0ph+nu0mh)*u(i1,i2,i3,vc) + nu0mh*u(i1-1,i2,i3,vc))*dxvsqi(0)\
112 ! +2.*(nu1ph*u(i1,i2+1,i3,vc) -(nu1ph+nu1mh)*u(i1,i2,i3,vc) + nu1mh*u(i1,i2-1,i3,vc))*dxvsqi(1)\
113 ! + ( nupzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))\
114 ! -numzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dx(0)*dx(1))
115 
116  if( computeRHS.eq.1 )then
117 
118  #If #dir == "0"
119  ! remove dir=0 implicit terms:
120  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu) + uu(uc)*dtScale/dt(i1,i2,i3) \
121  -u(i1,i2,i3,vc)*uy2(uc)-ux2(pc)\
122  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
123  \
124  + (nu1ph*u(i1,i2+1,i3,uc) + nu1mh*u(i1,i2-1,i3,uc))*dxvsqi(1)\
125  + (nuzpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))\
126  -nuzmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dx(0)*dx(1))
127 
128  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv) + uu(vc)*dtScale/dt(i1,i2,i3) \
129  -u(i1,i2,i3,vc)*uy2(vc)-uy2(pc)\
130  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
131  \
132  +2.*(nu1ph*u(i1,i2+1,i3,vc) + nu1mh*u(i1,i2-1,i3,vc))*dxvsqi(1)\
133  + (nupzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))\
134  -numzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dx(0)*dx(1))
135 
136  #Else
137 
138  ! remove dir=1 implicit terms:
139  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu) + uu(uc)*dtScale/dt(i1,i2,i3) \
140  -u(i1,i2,i3,uc)*ux2(uc) -ux2(pc)\
141  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
142  +2.*(nu0ph*u(i1+1,i2,i3,uc) + nu0mh*u(i1-1,i2,i3,uc))*dxvsqi(0)\
143  \
144  + (nuzpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))\
145  -nuzmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dx(0)*dx(1))
146 
147  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv) + uu(vc)*dtScale/dt(i1,i2,i3) \
148  -u(i1,i2,i3,uc)*ux2(vc) -uy2(pc)\
149  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
150  + (nu0ph*u(i1+1,i2,i3,vc) + nu0mh*u(i1-1,i2,i3,vc))*dxvsqi(0)\
151  \
152  + (nupzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))\
153  -numzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dx(0)*dx(1))
154 
155 
156  #End
157  if( use2ndOrderAD.eq.1 )then
158  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+ adE ## dir(i1,i2,i3,uc)
159  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+ adE ## dir(i1,i2,i3,vc)
160  end if
161  end if
162  else
163  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
164  am(i1,i2,i3)=0.
165  bm(i1,i2,i3)=1.
166  cm(i1,i2,i3)=0.
167  end if
168  if( computeRHS.eq.1 )then
169  f(i1,i2,i3,fcu)=uu(uc)
170  f(i1,i2,i3,fcv)=uu(vc)
171  end if
172  end if
173  endLoops()
174 
175 else if( nd.eq.3 )then
176 
177  stop 2945
178 
179 else
180  stop 888 ! unexpected value for nd
181 end if
182 
183 #endMacro
184 
185 
186 c ===================================================================================
187 c VP: incompressible Navier Stokes with a visco-plastic model
188 c *** curvilinear grid version ***
189 c Macro arguments:
190 c dir: 0,1,2
191 c====================================================================================
192 #beginMacro fillEquationsCurvilinearGridVP(dir)
193 
194  if( use4thOrderAD.eq.1 )then
195  write(*,*) 'insLineSolve: 4th order diss not finished'
196  stop 7654
197  end if
198 
199  ! write(*,'(" -- fill tridiagonal matrix curvlinear VP --")')
200 
201  ! set default values for no 2nd order artificial diffusion:
202  cdm=0.
203  cdDiag=0.
204  cdp=0.
205 
206 if( nd.eq.2 )then
207 
208  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
209  defineDerivativeMacros(2,2,curvilinear)
210 
211 
212  beginLoops()
213  if( mask(i1,i2,i3).gt.0 )then
214  if( use2ndOrderAD.eq.1 )then
215  defineArtificialDiffusionCoefficients(2,dir,C,)
216  end if
217 
218  ! Get the nonlinear viscosity at nearby points:
219  nuzmz=u(i1 ,i2-1,i3,vsc)
220  numzz=u(i1-1,i2 ,i3,vsc)
221  nuzzz=u(i1 ,i2 ,i3,vsc)
222  nupzz=u(i1+1,i2 ,i3,vsc)
223  nuzpz=u(i1 ,i2+1,i3,vsc)
224  ! Evaluate the nonlinear viscosity
225  ! getViscoPlasticViscosityCoefficient(nuzmz,i1 ,i2-1,i3,2,curvilinear)
226  ! getViscoPlasticViscosityCoefficient(numzz,i1-1,i2 ,i3,2,curvilinear)
227  ! getViscoPlasticViscosityCoefficient(nuzzz,i1 ,i2 ,i3,2,curvilinear)
228  ! getViscoPlasticViscosityCoefficient(nupzz,i1+1,i2 ,i3,2,curvilinear)
229  ! getViscoPlasticViscosityCoefficient(nuzpz,i1 ,i2+1,i3,2,curvilinear)
230 
231  ! u.t + u.grad(u) + p.x = Dx( 2*nu*u.x ) + Dy( nu*u.y ) + Dy( nu*v.x )
232  ! v.t + u.grad(v) + p.y = Dx( nu*v.x ) + Dy( 2*nu*v.y ) + Dx( nu*u.y )
233 
234  ! evaluate the jacobian at nearby points:
235  ajzmz = ajac2d(i1 ,i2-1,i3)
236  ajmzz = ajac2d(i1-1,i2 ,i3)
237  ajzzz = ajac2d(i1 ,i2 ,i3)
238  ajpzz = ajac2d(i1+1,i2 ,i3)
239  ajzpz = ajac2d(i1 ,i2+1,i3)
240 
241  ! 1. Get coefficients a11ph, a11mh, a22ph, etc. for
242  ! Dx( 2*nu*u.x ) + Dy( nu*u.y )
243  getCoeffForDxADxPlusDyBDy(au, 2.*nuzmz,2.*numzz,2.*nuzzz,2.*nupzz,2.*nuzpz, nuzmz,numzz,nuzzz,nupzz,nuzpz )
244 
245  ! 1b. Get coefficients b11ph,b11mh, etc. for
246  ! Dy( nu*v.x )
247  getCoeffForDyADx( bu, nuzmz,numzz,nuzzz,nupzz,nuzpz )
248 
249 
250  ! 2. Dx( nu*v.x ) + Dy( 2*nu*v.y )
251  getCoeffForDxADxPlusDyBDy(av, nuzmz,numzz,nuzzz,nupzz,nuzpz, 2.*nuzmz,2.*numzz,2.*nuzzz,2.*nupzz,2.*nuzpz )
252 
253  ! 2b. Dx( nu*u.y )
254  getCoeffForDxADy( bv, nuzmz,numzz,nuzzz,nupzz,nuzpz )
255 
256 
257  if( computeMatrix.eq.1 )then
258 
259  t1=(uu(uc)*rxi(dir,0)+uu(vc)*rxi(dir,1))*drv2i(dir)
260 
261  dr0i = 1./(ajzzz*dr(0)**2)
262  dr1i = 1./(ajzzz*dr(1)**2)
263 
264  #If #dir == "0"
265  if( systemComponent.eq.uc )then ! solve for u
266 
267  am(i1,i2,i3)= -t1 -cdm - au11mh*dr0i
268  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (au11ph+au11mh)*dr0i +(au22ph+au22mh)*dr1i
269  cm(i1,i2,i3)= t1 -cdp - au11ph*dr0i
270  else if( systemComponent.eq.vc )then ! solve for v
271  am(i1,i2,i3)= -t1 -cdm - av11mh*dr0i
272  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (av11ph+av11mh)*dr0i +(av22ph+av22mh)*dr1i
273  cm(i1,i2,i3)= t1 -cdp - av11ph*dr0i
274  else
275  stop 6139
276  end if
277  #Else
278  if( systemComponent.eq.uc )then ! solve for u
279  am(i1,i2,i3)= -t1 -cdm - au22mh*dr1i
280  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (au11ph+au11mh)*dr0i +(au22ph+au22mh)*dr1i
281  cm(i1,i2,i3)= t1 -cdp - au22ph*dr1i
282  else if( systemComponent.eq.vc )then ! solve for v
283  am(i1,i2,i3)= -t1 -cdm - av22mh*dr1i
284  bm(i1,i2,i3)= dtScale/dt(i1,i2,i3) +cdDiag + (av11ph+av11mh)*dr0i +(av22ph+av22mh)*dr1i
285  cm(i1,i2,i3)= t1 -cdp - av22ph*dr1i
286  else
287  stop 6139
288  end if
289  #End
290 
291  end if
292 
293 ! Here is the full residual:
294 !
295 ! 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)\
296 ! -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
297 ! + ( \
298 ! + ( au11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - au11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
299 ! + ( au22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - au22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
300 ! + (au12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-au12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
301 ! + (au21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-au21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
302 ! \
303 ! + ( bu11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - bu11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
304 ! + ( bu22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - bu22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
305 ! + (bu12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-bu12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
306 ! + (bu21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-bu21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
307 ! )/ajzzz
308 !
309 ! 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)\
310 ! -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
311 ! + ( \
312 ! + ( av11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - av11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
313 ! + ( av22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - av22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
314 ! + (av12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-av12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
315 ! + (av21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-av21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
316 ! \
317 ! + ( bv11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - bv11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
318 ! + ( bv22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - bv22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
319 ! + (bv12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-bv12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
320 ! + (bv21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-bv21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
321 ! )/ajzzz
322 
323 
324 
325  if( computeRHS.eq.1 )then
326 
327  #If #dir == "0"
328  ! remove dir=0 implicit terms:
329  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu) + uu(uc)*dtScale/dt(i1,i2,i3) \
330  +(-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1))*us(uc) -ux2c(pc)\
331  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
332  + ( \
333  \
334  + ( au22ph*(u(i1,i2+1,i3,uc) ) - au22mh*( -u(i1,i2-1,i3,uc)) )/dr(1)**2\
335  + (au12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-au12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
336  + (au21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-au21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
337  \
338  + ( bu11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - bu11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
339  + ( bu22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - bu22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
340  + (bu12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-bu12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
341  + (bu21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-bu21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
342  )/ajzzz
343 
344  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv) + uu(vc)*dtScale/dt(i1,i2,i3) \
345  +(-uu(uc)*rxi(1,0)-uu(vc)*rxi(1,1))*us(vc) -uy2c(pc)\
346  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
347  + ( \
348  \
349  + ( av22ph*(u(i1,i2+1,i3,vc) ) - av22mh*( -u(i1,i2-1,i3,vc)) )/dr(1)**2\
350  + (av12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-av12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
351  + (av21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-av21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
352  \
353  + ( bv11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - bv11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
354  + ( bv22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - bv22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
355  + (bv12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-bv12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
356  + (bv21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-bv21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
357  )/ajzzz
358 
359  #Else
360 
361  ! remove dir=1 implicit terms:
362  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu) + uu(uc)*dtScale/dt(i1,i2,i3) \
363  + (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1))*ur(uc) -ux2c(pc) \
364  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
365  + ( \
366  + ( au11ph*(u(i1+1,i2,i3,uc) ) - au11mh*( -u(i1-1,i2,i3,uc)) )/dr(0)**2\
367  \
368  + (au12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-au12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
369  + (au21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-au21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
370  \
371  + ( bu11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - bu11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
372  + ( bu22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - bu22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
373  + (bu12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-bu12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
374  + (bu21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-bu21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
375  )/ajzzz
376 
377  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv) + uu(vc)*dtScale/dt(i1,i2,i3) \
378  + (-uu(uc)*rxi(0,0)-uu(vc)*rxi(0,1))*ur(vc) -uy2c(pc)\
379  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
380  + ( \
381  + ( av11ph*(u(i1+1,i2,i3,vc) ) - av11mh*( -u(i1-1,i2,i3,vc)) )/dr(0)**2\
382  \
383  + (av12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-av12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
384  + (av21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-av21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
385  \
386  + ( bv11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - bv11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
387  + ( bv22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - bv22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
388  + (bv12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-bv12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
389  + (bv21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-bv21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
390  )/ajzzz
391 
392  #End
393  if( use2ndOrderAD.eq.1 )then
394  f(i1,i2,i3,fcu)=f(i1,i2,i3,fcu)+ adE ## dir(i1,i2,i3,uc)
395  f(i1,i2,i3,fcv)=f(i1,i2,i3,fcv)+ adE ## dir(i1,i2,i3,vc)
396  end if
397  end if
398  else
399  if( computeMatrix.eq.1 )then ! for interpolation points or unused:
400  am(i1,i2,i3)=0.
401  bm(i1,i2,i3)=1.
402  cm(i1,i2,i3)=0.
403  end if
404  if( computeRHS.eq.1 )then
405  f(i1,i2,i3,fcu)=uu(uc)
406  f(i1,i2,i3,fcv)=uu(vc)
407  end if
408  end if
409  endLoops()
410 
411 else if( nd.eq.3 )then
412 
413  stop 2945
414 
415 else
416  stop 888 ! unexpected value for nd
417 end if
418 
419 #endMacro
420 
421 
422 
423 
424 c==================================================================================
425 c Residual Computation for VP: incompressible Navier Stokes with a visco-plastic model
426 c
427 c Macro args:
428 c GRIDTYPE: rectangular, curvilinear
429 c====================================================================================
430 #beginMacro computeResidualVP(GRIDTYPE)
431 
432  ! write(*,'("new: computeResidualVP, use2ndOrderAD=",i2)') use2ndOrderAD
433 
434  if( use4thOrderAD.eq.1 )then
435  write(*,*) 'insLineSolve: computeResidualVP: 4th order diss not finished'
436  stop 7654
437  end if
438 
439  if( nd.eq.2 )then
440 
441  ! defineDerivativeMacros(DIM,ORDER,GRIDTYPE) : defineMacro UX(cc) ux22r(i1,i2,i3,cc) etc.
442  ! defineDerivativeMacros(2,2,GRIDTYPE)
443 
444  beginLoops()
445  if( mask(i1,i2,i3).gt.0 )then
446 
447  ! Get the nonlinear viscosity at nearby points:
448  nuzmz=u(i1 ,i2-1,i3,vsc)
449  numzz=u(i1-1,i2 ,i3,vsc)
450  nuzzz=u(i1 ,i2 ,i3,vsc)
451  nupzz=u(i1+1,i2 ,i3,vsc)
452  nuzpz=u(i1 ,i2+1,i3,vsc)
453  ! Evaluate the nonlinear viscosity "nu"
454  ! getViscoPlasticViscosityCoefficient(nuzmz,i1 ,i2-1,i3,2,GRIDTYPE)
455  ! getViscoPlasticViscosityCoefficient(numzz,i1-1,i2 ,i3,2,GRIDTYPE)
456  ! getViscoPlasticViscosityCoefficient(nuzzz,i1 ,i2 ,i3,2,GRIDTYPE)
457  ! getViscoPlasticViscosityCoefficient(nupzz,i1+1,i2 ,i3,2,GRIDTYPE)
458  ! getViscoPlasticViscosityCoefficient(nuzpz,i1 ,i2+1,i3,2,GRIDTYPE)
459 
460  #If #GRIDTYPE == "rectangular"
461  nu0ph = .5*( nupzz+nuzzz ) ! nu(i1+1/2,i2,i3)
462  nu0mh = .5*( nuzzz+numzz ) ! nu(i1-1/2,i2,i3)
463 
464  nu1ph = .5*( nuzpz+nuzzz ) ! nu(i1,i2+1/2,i3)
465  nu1mh = .5*( nuzzz+nuzmz ) ! nu(i1,i2-1/2,i3)
466 
467  ! u.t + u.grad(u) + p.x = Dx( 2*nu*u.x ) + Dy( nu*u.y ) + Dy( nu*v.x )
468  ! v.t + u.grad(v) + p.y = Dx( nu*v.x ) + Dy( 2*nu*v.y ) + Dx( nu*u.y )
469  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)\
470  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
471  +2.*(nu0ph*u(i1+1,i2,i3,uc) -(nu0ph+nu0mh)*u(i1,i2,i3,uc) + nu0mh*u(i1-1,i2,i3,uc))*dxvsqi(0)\
472  + (nu1ph*u(i1,i2+1,i3,uc) -(nu1ph+nu1mh)*u(i1,i2,i3,uc) + nu1mh*u(i1,i2-1,i3,uc))*dxvsqi(1)\
473  + (nuzpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))\
474  -nuzmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dx(0)*dx(1))
475 
476  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)\
477  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
478  + (nu0ph*u(i1+1,i2,i3,vc) -(nu0ph+nu0mh)*u(i1,i2,i3,vc) + nu0mh*u(i1-1,i2,i3,vc))*dxvsqi(0)\
479  +2.*(nu1ph*u(i1,i2+1,i3,vc) -(nu1ph+nu1mh)*u(i1,i2,i3,vc) + nu1mh*u(i1,i2-1,i3,vc))*dxvsqi(1)\
480  + (nupzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))\
481  -numzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dx(0)*dx(1))
482 
483  if( use2ndOrderAD.eq.1 )then
484  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint2dR(i1,i2,i3,uc)
485  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint2dR(i1,i2,i3,vc)
486  end if
487 
488  if( computeTemperature.ne.0 )then
489  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)
490  ! --- artificial dissipation for T: do this for now:
491  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
492  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint2dR(i1,i2,i3,tc)
493  else if( use4thOrderAD.eq.1 )then
494  ! compute adc2, adc4:
495  getDerivativesAndDissipation(2,2,GRIDTYPE)
496  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad2(adc2,tc)+ad4(adc4,tc)
497  end if
498  end if
499 
500  #Else
501  ! ************ VP curvilinear case ********************
502 
503  ! evaluate the jacobian at nearby points:
504  ajzmz = ajac2d(i1 ,i2-1,i3)
505  ajmzz = ajac2d(i1-1,i2 ,i3)
506  ajzzz = ajac2d(i1 ,i2 ,i3)
507  ajpzz = ajac2d(i1+1,i2 ,i3)
508  ajzpz = ajac2d(i1 ,i2+1,i3)
509 
510  ! 1. Get coefficients a11ph, a11mh, a22ph, etc. for
511  ! Dx( 2*nu*u.x ) + Dy( nu*u.y )
512  getCoeffForDxADxPlusDyBDy(a, 2.*nuzmz,2.*numzz,2.*nuzzz,2.*nupzz,2.*nuzpz, nuzmz,numzz,nuzzz,nupzz,nuzpz )
513 
514  ! 1b. Get coefficients b11ph,b11mh, etc. for
515  ! Dy( nu*v.x )
516  getCoeffForDyADx( b, nuzmz,numzz,nuzzz,nupzz,nuzpz )
517 
518 
519  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)\
520  -thermalExpansivity*gravity(0)*u(i1,i2,i3,tc) \
521  + ( \
522  + ( a11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - a11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
523  + ( a22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - a22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
524  + (a12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-a12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
525  + (a21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-a21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
526  \
527  + ( b11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - b11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
528  + ( b22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - b22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
529  + (b12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-b12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
530  + (b21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-b21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
531  )/ajzzz
532 
533  ! 2. Dx( nu*v.x ) + Dy( 2*nu*v.y )
534  getCoeffForDxADxPlusDyBDy(a, nuzmz,numzz,nuzzz,nupzz,nuzpz, 2.*nuzmz,2.*numzz,2.*nuzzz,2.*nupzz,2.*nuzpz )
535 
536  ! 2b. Dx( nu*u.y )
537  getCoeffForDxADy( b, nuzmz,numzz,nuzzz,nupzz,nuzpz )
538 
539 
540  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)\
541  -thermalExpansivity*gravity(1)*u(i1,i2,i3,tc) \
542  + ( \
543  + ( a11ph*(u(i1+1,i2,i3,vc)-u(i1,i2,i3,vc)) - a11mh*(u(i1,i2,i3,vc)-u(i1-1,i2,i3,vc)) )/dr(0)**2\
544  + ( a22ph*(u(i1,i2+1,i3,vc)-u(i1,i2,i3,vc)) - a22mh*(u(i1,i2,i3,vc)-u(i1,i2-1,i3,vc)) )/dr(1)**2\
545  + (a12pzz*(u(i1+1,i2+1,i3,vc)-u(i1+1,i2-1,i3,vc))-a12mzz*(u(i1-1,i2+1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
546  + (a21zpz*(u(i1+1,i2+1,i3,vc)-u(i1-1,i2+1,i3,vc))-a21zmz*(u(i1+1,i2-1,i3,vc)-u(i1-1,i2-1,i3,vc)))/(4.*dr(0)*dr(1))\
547  \
548  + ( b11ph*(u(i1+1,i2,i3,uc)-u(i1,i2,i3,uc)) - b11mh*(u(i1,i2,i3,uc)-u(i1-1,i2,i3,uc)) )/dr(0)**2\
549  + ( b22ph*(u(i1,i2+1,i3,uc)-u(i1,i2,i3,uc)) - b22mh*(u(i1,i2,i3,uc)-u(i1,i2-1,i3,uc)) )/dr(1)**2\
550  + (b12pzz*(u(i1+1,i2+1,i3,uc)-u(i1+1,i2-1,i3,uc))-b12mzz*(u(i1-1,i2+1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
551  + (b21zpz*(u(i1+1,i2+1,i3,uc)-u(i1-1,i2+1,i3,uc))-b21zmz*(u(i1+1,i2-1,i3,uc)-u(i1-1,i2-1,i3,uc)))/(4.*dr(0)*dr(1))\
552  )/ajzzz
553 
554  if( use2ndOrderAD.eq.1 )then
555  residual(i1,i2,i3,uc)=residual(i1,i2,i3,uc)+adSelfAdjoint2dC(i1,i2,i3,uc)
556  residual(i1,i2,i3,vc)=residual(i1,i2,i3,vc)+adSelfAdjoint2dC(i1,i2,i3,vc)
557  end if
558 
559  if( computeTemperature.ne.0 )then
560  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)
561  ! --- artificial dissipation for T: do this for now:
562  if( use2ndOrderAD.eq.1 .and. use4thOrderAD.eq.0 )then
563  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+adSelfAdjoint2dC(i1,i2,i3,tc)
564  else if( use4thOrderAD.eq.1 )then
565  ! compute adc2, adc4:
566  getDerivativesAndDissipation(2,2,GRIDTYPE)
567  residual(i1,i2,i3,tc)=residual(i1,i2,i3,tc)+ad2(adc2,tc)+ad4(adc4,tc)
568  end if
569  end if
570 
571  #End
572 
573  end if
574  endLoops()
575 
576  else if( nd.eq.3 )then
577 
578  stop 2945
579 
580  else
581  stop 888 ! unexpected value for nd
582  end if
583 
584 #endMacro
Generated on Fri Jan 4 2013 10:19:24 for CG by   doxygen 1.8.3