Overture  Version 25
neumannEquationForcing.h
Go to the documentation of this file.
1 ! This Ogmg macro defines the forcing for fourth-order Neumann boundary conditions
2 
3 ! define derivatives in the r, s, and t directions:
4 
5 ! 2nd-order centered:
6 #defineMacro FR(i1,i2,i3) ((f(i1+1,i2,i3)-f(i1-1,i2,i3))*D12(0))
7 #defineMacro FS(i1,i2,i3) ((f(i1,i2+1,i3)-f(i1,i2-1,i3))*D12(1))
8 #defineMacro FT(i1,i2,i3) ((f(i1,i2,i3+1)-f(i1,i2,i3-1))*D12(2))
9 
10 #defineMacro FRR(i1,i2,i3) ((f(i1+1,i2,i3)-2.*f(i1,i2,i3)+f(i1-1,i2,i3))*D22(0))
11 #defineMacro FSS(i1,i2,i3) ((f(i1,i2+1,i3)-2.*f(i1,i2,i3)+f(i1,i2-1,i3))*D22(1))
12 #defineMacro FTT(i1,i2,i3) ((f(i1,i2,i3+1)-2.*f(i1,i2,i3)+f(i1,i2,i3-1))*D22(2))
13 
14 
15 ! one-sided approximations on the left:
16 #defineMacro FRa(i1,i2,i3) ((-f(i1+2,i2,i3)+4.*f(i1+1,i2,i3)-3.*f(i1,i2,i3))*D12(0))
17 #defineMacro FSa(i1,i2,i3) ((-f(i1,i2+2,i3)+4.*f(i1,i2+1,i3)-3.*f(i1,i2,i3))*D12(1))
18 #defineMacro FTa(i1,i2,i3) ((-f(i1,i2,i3+2)+4.*f(i1,i2,i3+1)-3.*f(i1,i2,i3))*D12(2))
19 
20 #defineMacro FRRa(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1+1,i2,i3)+4.*f(i1+2,i2,i3)-f(i1+3,i2,i3))*D22(0))
21 #defineMacro FSSa(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1,i2+1,i3)+4.*f(i1,i2+2,i3)-f(i1,i2+3,i3))*D22(1))
22 #defineMacro FTTa(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1,i2,i3+1)+4.*f(i1,i2,i3+2)-f(i1,i2,i3+3))*D22(2))
23 
24 ! one-sided approximations on the right:
25 #defineMacro FRb(i1,i2,i3) (( f(i1-2,i2,i3)-4.*f(i1-1,i2,i3)+3.*f(i1,i2,i3))*D12(0))
26 #defineMacro FSb(i1,i2,i3) (( f(i1,i2-2,i3)-4.*f(i1,i2-1,i3)+3.*f(i1,i2,i3))*D12(1))
27 #defineMacro FTb(i1,i2,i3) (( f(i1,i2,i3-2)-4.*f(i1,i2,i3-1)+3.*f(i1,i2,i3))*D12(2))
28 
29 #defineMacro FRRb(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1-1,i2,i3)+4.*f(i1-2,i2,i3)-f(i1-3,i2,i3))*D22(0))
30 #defineMacro FSSb(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1,i2-1,i3)+4.*f(i1,i2-2,i3)-f(i1,i2-3,i3))*D22(1))
31 #defineMacro FTTb(i1,i2,i3) ((2.*f(i1,i2,i3)-5.*f(i1,i2,i3-1)+4.*f(i1,i2,i3-2)-f(i1,i2,i3-3))*D22(2))
32 
33 ! 2nd-order centered:
34 #defineMacro GVR() ((gv(+1,0,0)-gv(-1,0,0))*D12(0))
35 #defineMacro GVS() ((gv(0,+1,0)-gv(0,-1,0))*D12(1))
36 #defineMacro GVT() ((gv(0,0,+1)-gv(0,0,-1))*D12(2))
37 
38 #defineMacro GVRR() ((gv(+1,0,0)-2.*gv(0,0,0)+gv(-1,0,0))*D22(0))
39 #defineMacro GVSS() ((gv(0,+1,0)-2.*gv(0,0,0)+gv(0,-1,0))*D22(1))
40 #defineMacro GVTT() ((gv(0,0,+1)-2.*gv(0,0,0)+gv(0,0,-1))*D22(2))
41 
42 #defineMacro GVRS() ( ((gv(+1,+1,0)-gv(+1,-1,0)) - (gv(-1,+1,0)-gv(-1,-1,0)))*d12(0)*d12(1) )
43 #defineMacro GVRT() ( ((gv(+1,0,+1)-gv(+1,0,-1)) - (gv(-1,0,+1)-gv(-1,0,-1)))*d12(0)*d12(2) )
44 #defineMacro GVST() ( ((gv(0,+1,+1)-gv(0,+1,-1)) - (gv(0,-1,+1)-gv(0,-1,-1)))*d12(1)*d12(2) )
45 
46 
47 #defineMacro FVR() ((fv(+1,0,0)-fv(-1,0,0))*D12(0))
48 #defineMacro FVS() ((fv(0,+1,0)-fv(0,-1,0))*D12(1))
49 #defineMacro FVT() ((fv(0,0,+1)-fv(0,0,-1))*D12(2))
50 
51 #defineMacro FVRR() ((fv(+1,0,0)-2.*fv(0,0,0)+fv(-1,0,0))*D22(0))
52 #defineMacro FVSS() ((fv(0,+1,0)-2.*fv(0,0,0)+fv(0,-1,0))*D22(1))
53 #defineMacro FVTT() ((fv(0,0,+1)-2.*fv(0,0,0)+fv(0,0,-1))*D22(2))
54 
55 ! Extrapolation of f(i1,i2,i3) in direction (is1,is2,is3)
56 #defineMacro extrap1(f,i1,i2,i3,is1,is2,is3) (f(i1+(is1),i2+(is2),i3+(is3)))
57 
58 #defineMacro extrap2(f,i1,i2,i3,is1,is2,is3) (2.*f(i1+(is1),i2+(is2),i3+(is3))-f(i1+2*(is1),i2+2*(is2),i3+2*(is3)))
59 
60 #defineMacro extrap3(f,i1,i2,i3,is1,is2,is3) (3.*f(i1+(is1),i2+(is2),i3+(is3))-3.*f(i1+2*(is1),i2+2*(is2),i3+2*(is3))+f(i1+3*(is1),i2+3*(is2),i3+3*(is3)))
61 
62 #defineMacro extrap4(f,i1,i2,i3,is1,is2,is3) (4.*f(i1+(is1),i2+(is2),i3+(is3))-6.*f(i1+2*(is1),i2+2*(is2),i3+2*(is3))+4.*f(i1+3*(is1),i2+3*(is2),i3+3*(is3))-f(i1+4*(is1),i2+4*(is2),i3+4*(is3)))
63 
64 
65 ! ================================================================================================
66 ! Extrapolate a point. Choose the order of extrapolation based on how many
67 ! valid points exist (mask>0)
68 !
69 ! fe : put result here
70 ! (k1,k2,k3) : check mask at points (i1+m*is1,i1+m*is2,i3+m*is3) m=1,2,..
71 ! (l1,l2,l3) : Extrapolate point (l1,l2,l3) using points (l1+m*is1,l1+m*is2,l3+m*is3)
72 ! (is1,is2,is3) : direction (shift) of extrapolation
73 !
74 ! ================================================================================================
75 #beginMacro extrapWithMask(fe, f, k1,k2,k3, l1,l2,l3, is1,is2,is3 )
76  if( mask(k1+ (is1),k2+ (is2),k3+ (is3)).gt.0 .and. \
77  mask(k1+2*(is1),k2+2*(is2),k3+2*(is3)).gt.0 .and.\ mask(k1+3*(is1),k2+3*(is2),k3+3*(is3)).gt.0 .and.\ mask(k1+4*(is1),k2+4*(is2),k3+4*(is3)).gt.0 )then
78  fe=extrap4(f,l1,l2,l3,is1,is2,is3)
79  else if( mask(k1+ (is1),k2+ (is2),k3+ (is3)).gt.0 .and. \
80  mask(k1+2*(is1),k2+2*(is2),k3+2*(is3)).gt.0 .and.\ mask(k1+3*(is1),k2+3*(is2),k3+3*(is3)).gt.0 )then
81  fe=extrap3(f,l1,l2,l3,is1,is2,is3)
82  else if( mask(k1+ (is1),k2+ (is2),k3+ (is3)).gt.0 .and. \
83  mask(k1+2*(is1),k2+2*(is2),k3+2*(is3)).gt.0 )then
84  fe=extrap2(f,l1,l2,l3,is1,is2,is3)
85  else
86  fe=extrap1(f,l1,l2,l3,is1,is2,is3)
87  end if
88 #endMacro
89 
90 ! ===================================================================================================
91 ! Macro: Evaluate the boundary forcing g at the ghost point (k1,k2,k3) next to the target
92 ! ghost point (l1,l2,l3) (where we are evaluating derivatives of g) and the boundary point (b1,b2,b3)
93 !
94 ! | | |
95 ! ---B---L---X---- <- boundary
96 ! | | |
97 ! ---K---+---+---- <- ghost
98 !
99 ! This macro expects the following to be already set:
100 ! ax1, ax2 : tangential directions:
101 ! ax1= mod(axis+1,nd)
102 ! ax2= mod(axis+2,nd)
103 ! mdim(0:1,0:2) : index bounds on valid points where g is defined.
104 ! ===================================================================================================
105 #beginMacro getAdjacentBoundaryForcing(k1,k2,k3, l1,l2,l3, b1,b2,b3)
106 
107  ! Add these checks -- comment out later
108  if( abs(k1-l1).gt.1 .or. abs(k2-l2).gt.1 .or. abs(k3-l3).gt.1 )then
109  stop 3338
110  end if
111  ! iv(0:2) = boundary point next to (k1,k2,k3)
112  iv(0)=b1
113  iv(1)=b2
114  iv(2)=b3
115  ! dv(0:2) : extrapolate in direction dv (if dv[a]=0 a=0,1,2 then there is no extrapolation)
116  dv(0)=0
117  dv(1)=0
118  dv(2)=0
119  if( iv(ax1).lt.mdim(0,ax1) )then
120  dv(ax1)=1
121  else if( iv(ax1).gt.mdim(1,ax1) )then
122  dv(ax1)=-1
123  end if
124  if( iv(ax2).lt.mdim(0,ax2) )then
125  dv(ax2)=1
126  else if( iv(ax2).gt.mdim(1,ax2) )then
127  dv(ax2)=-1
128  end if
129  if( mask(iv(0)+dv(0),iv(1)+dv(1),iv(2)+dv(2)).le.0 )then
130  ! The neighbouring pt kv or the point we start the extrapolation from is not valid
131  ! Find a direction to extrapolate in: The target point must be valid so start the
132  ! extrpolation in that direction.
133  dv(0)=l1-k1
134  dv(1)=l2-k2
135  dv(2)=l3-k3
136  end if
137  if( dv(0).eq.0 .and. dv(1).eq.0. .and. dv(2).eq.0 )then
138  ! We can use the value from the adjacent point:
139  gv(k1-l1,k2-l2,k3-l3)=f(k1,k2,k3)
140  else if( mask(iv(0)+ dv(0),iv(1)+ dv(1),iv(2)+ dv(2)).gt.0 .and. \
141  mask(iv(0)+2*dv(0),iv(1)+2*dv(1),iv(2)+2*dv(2)).gt.0 .and. \
142  mask(iv(0)+3*dv(0),iv(1)+3*dv(1),iv(2)+3*dv(2)).gt.0 .and. \
143  mask(iv(0)+4*dv(0),iv(1)+4*dv(1),iv(2)+4*dv(2)).gt.0 )then
144  ! we can extrapolate with 4 points:
145  gv(k1-l1,k2-l2,k3-l3) = extrap4(f,k1,k2,k3, dv(0), dv(1), dv(2))
146  else if( mask(iv(0)+ dv(0),iv(1)+ dv(1),iv(2)+ dv(2)).gt.0 .and. \
147  mask(iv(0)+2*dv(0),iv(1)+2*dv(1),iv(2)+2*dv(2)).gt.0 .and. \
148  mask(iv(0)+3*dv(0),iv(1)+3*dv(1),iv(2)+3*dv(2)).gt.0 )then
149  ! we can extrapolate with 3 points:
150  gv(k1-l1,k2-l2,k3-l3) = extrap3(f,k1,k2,k3, dv(0), dv(1), dv(2))
151  else if( mask(iv(0)+ dv(0),iv(1)+ dv(1),iv(2)+ dv(2)).gt.0 .and. \
152  mask(iv(0)+2*dv(0),iv(1)+2*dv(1),iv(2)+2*dv(2)).gt.0 )then
153  ! we can extrapolate with 2 points:
154  gv(k1-l1,k2-l2,k3-l3) = extrap2(f,k1,k2,k3, dv(0), dv(1), dv(2))
155  else
156  ! as a backup just use the value from the target point
157  gv(k1-l1,k2-l2,k3-l3)=f(l1,l2,l3)
158  end if
159 #endMacro
160 
161 
162 ! ========================================================================================================================
163 ! This Ogmg macro defines the forcing for fourth-order Neumann boundary conditions
164 ! Lu = ff , Bu=g
165 ! Input:
166 ! [mm1a,mm1b][mm2a,mm2b][mm3a,mm3b] : indexes for the boundary
167 ! (i1,i2,i3) : point on the boundary
168 ! (j1,j2,j3) : ghost point
169 ! (is1,is2,is3) : usual
170 ! FORCING: forcing or no forcing (if noForcing then return values are all zero)
171 ! GRIDTYPE = rectangular or curvilinear
172 ! DIR = R or S or T
173 ! DIM = 2 or 3
174 ! Return: DIM=2, DIR=R : ff, ffr, g, gss
175 ! DIM=2, DIR=S : ff, ffs, g, grr
176 ! DIM=3, DIR=R : ff, ffr, g, gss, gtt + curvilinear: ffs, fft, gst
177 ! DIM=3, DIR=S : ff, ffs, g, grr, gtt + curvilinear: ffr, fft, grt
178 ! DIM=3, DIR=T : ff, fft, g, grr, gss + curvilinear: ffr, ffs, grs
179 ! ========================================================================================================================
180 #beginMacro defineNeumannEquationForcing(mm1a,mm1b,mm2a,mm2b,mm3a,mm3b,FORCING,GRIDTYPE,DIR,DIM)
181  ! the rhs for the mixed BC is stored in the ghost point value of f
182 #If #GRIDTYPE eq "rectangular"
183 
184 ! Cartesian grids use dx:
185 #defineMacro D12(axis) h12(axis)
186 #defineMacro D22(axis) h22(axis)
187 
188  #If #FORCING eq "forcing"
189  g = f(j1,j2,j3)
190  ff=f(i1,i2,i3)
191 
192  #If #DIR eq "R"
193  ! Note "g" is located on the ghost point "j1" of f
194 
195  ! 2nd-order one sided:
196  ! ffr=(-f(i1+2*is1,i2,i3)+4.*f(i1+is1,i2,i3)-3.*ff)/(2.*dx(0))
197  ! 3rd-order one sided: 100510 -- added is1
198  ffr=is1*(-11.*ff+18.*f(i1+is1,i2,i3)-9.*f(i1+2*is1,i2,i3)+2.*f(i1+3*is1,i2,i3))/(6.*dx(0))
199 
200  ! 100610: Check the mask for computing valid tangential derivatives:
201 
202  ! NOTE: the forcing f and g are only assumed to be given where mask>0
203  ! In order to compute tangential derivatives of the forcing we may need to fill in
204  ! neighbouring values of the forcing at interp and unused points
205  gv( 0, 0, 0)=f(j1,i2,i3)
206  i2m1 = i2-1
207  if( i2m1.lt.mm2a .or. mask(i1,i2m1,i3).le.0 )then
208  ! f(j1,i2m1,i3)= extrap3(f,j1,i2m1,i3, 0,1,0)
209  ! gv( 0,-1, 0)=extrap3(f,j1,i2m1,i3, 0,1,0)
210  ! extrapWithMask
211  extrapWithMask( gv( 0,-1, 0), f, i1,i2m1,i3, j1,i2m1,i3, 0,1,0 )
212  else
213  gv( 0,-1, 0)=f(j1,i2m1,i3)
214  end if
215  i2p1 = i2+1
216  if( i2p1.gt.mm2b .or. mask(i1,i2p1,i3).le.0 )then
217  ! f(j1,i2p1,i3)= extrap3(f,j1,i2p1,i3, 0,-1,0)
218  ! gv( 0,+1, 0)=extrap3(f,j1,i2p1,i3, 0,-1,0)
219  extrapWithMask(gv( 0,+1, 0), f, i1,i2p1,i3, j1,i2p1,i3, 0,-1,0 )
220  else
221  gv( 0,+1, 0)=f(j1,i2p1,i3)
222  end if
223  ! gss=FSS(j1,i2,i3)
224  gss = GVSS()
225 
226  #If #DIM eq "3"
227  i3m1 = i3-1
228  if( i3m1.lt.mm3a .or. mask(i1,i2,i3m1).le.0 )then
229  ! f(j1,i2,i3m1)= extrap3(f,j1,i2,i3m1, 0,0,1)
230  ! gv( 0, 0,-1) = extrap3(f,j1,i2,i3m1, 0,0,1)
231  extrapWithMask(gv( 0, 0,-1), f, i1,i2,i3m1, j1,i2,i3m1, 0,0,1 )
232  else
233  gv( 0, 0,-1) = f(j1,i2,i3m1)
234  end if
235  i3p1 = i3+1
236  if( i3p1.gt.mm3b .or. mask(i1,i2,i3p1).le.0 )then
237  ! f(j1,i2,i3p1)= extrap3(f,j1,i2,i3p1, 0,0,-1)
238  ! gv( 0, 0,+1) = extrap3(f,j1,i2,i3p1, 0,0,-1)
239  extrapWithMask(gv( 0, 0,+1), f, i1,i2,i3p1, j1,i2,i3p1, 0,0,-1 )
240  else
241  gv( 0, 0,+1) = f(j1,i2,i3p1)
242  end if
243  ! gtt=FTT(j1,i2,i3)
244  gtt = GVTT()
245 
246  #End
247 
248  #Elif #DIR eq "S"
249  ! 2nd-order one sided:
250  ! ffs=(-f(i1,i2+2*is2,i3)+4.*f(i1,i2+is2,i3)-3.*ff)/(2.*dx(1))
251  ! 3rd-order one sided:
252  ffs=is2*(-11.*ff+18.*f(i1,i2+is2,i3)-9.*f(i1,i2+2*is2,i3)+2.*f(i1,i2+3*is2,i3))/(6.*dx(1))
253 
254  ! NOTE: the forcing f and g are only assumed to be given where mask>0
255  ! In order to compute tangential derivatives of the forcing we may need to fill in
256  ! neighbouring values of the forcing at interp and unused points
257  gv( 0, 0, 0)=f(i1,j2,i3)
258  i1m1 = i1-1
259  if( i1m1.lt.mm1a .or. mask(i1m1,i2,i3).le.0 )then
260  ! f(i1m1,j2,i3)= extrap3(f,i1m1,j2,i3, 1,0,0)
261  ! gv(-1, 0, 0) = extrap3(f,i1m1,j2,i3, 1,0,0)
262  extrapWithMask(gv(-1, 0, 0), f, i1m1,i2,i3, i1m1,j2,i3, 1,0,0 )
263  else
264  gv(-1, 0, 0) = f(i1m1,j2,i3)
265  end if
266  i1p1 = i1+1
267  if( i1p1.gt.mm1b .or. mask(i1p1,i2,i3).le.0 )then
268  ! f(i1p1,j2,i3)= extrap3(f,i1p1,j2,i3,-1,0,0)
269  ! gv(+1, 0, 0) = extrap3(f,i1p1,j2,i3,-1,0,0)
270  extrapWithMask(gv(+1, 0, 0), f, i1p1,i2,i3, i1p1,j2,i3, -1,0,0 )
271  else
272  gv(+1, 0, 0) = f(i1p1,j2,i3)
273  end if
274  ! grr=FRR(i1,j2,i3)
275  grr = GVRR()
276 
277  #If #DIM eq "3"
278  i3m1 = i3-1
279  if( i3m1.lt.mm3a .or. mask(i1,i2,i3m1).le.0 )then
280  ! f(i1,j2,i3m1)= extrap3(f,i1,j2,i3m1, 0,0,1)
281  ! gv( 0, 0,-1) = extrap3(f,i1,j2,i3m1, 0,0,1)
282  extrapWithMask(gv( 0, 0,-1), f, i1,i2,i3m1, i1,j2,i3m1, 0,0,1 )
283  else
284  gv( 0, 0,-1) = f(i1,j2,i3m1)
285  end if
286  i3p1 = i3+1
287  if( i3p1.gt.mm3b .or. mask(i1,i2,i3p1).le.0 )then
288  ! f(i1,j2,i3p1)= extrap3(f,i1,j2,i3p1, 0,0,-1)
289  ! gv( 0, 0,+1) = extrap3(f,i1,j2,i3p1, 0,0,-1)
290  extrapWithMask(gv( 0, 0,+1), f, i1,i2,i3p1, i1,j2,i3p1, 0,0,-1 )
291  else
292  gv( 0, 0,+1) = f(i1,j2,i3p1)
293  end if
294  ! gtt=FTT(i1,j2,i3)
295  gtt = GVTT()
296 
297  #End
298 
299  #Elif #DIR eq "T"
300  ! 3rd-order one sided:
301  fft=is3*(-11.*ff+18.*f(i1,i2,i3+is3)-9.*f(i1,i2,i3+2*is3)+2.*f(i1,i2,i3+3*is3))/(6.*dx(2))
302 
303  gv( 0, 0, 0)=f(i1,i2,j3)
304  i1m1 = i1-1
305  if( i1m1.lt.mm1a .or. mask(i1m1,i2,i3).le.0 )then
306  ! f(i1m1,i2,j3)= extrap3(f,i1m1,i2,j3, 1,0,0)
307  ! gv(-1, 0, 0) = extrap3(f,i1m1,i2,j3, 1,0,0)
308  extrapWithMask(gv(-1, 0, 0), f, i1m1,i2,i3, i1m1,i2,j3, 1,0,0 )
309  else
310  gv(-1, 0, 0) = f(i1m1,i2,j3)
311  end if
312  i1p1 = i1+1
313  if( i1p1.gt.mm1b .or. mask(i1p1,i2,i3).le.0 )then
314  ! f(i1p1,i2,j3)= extrap3(f,i1p1,i2,j3,-1,0,0)
315  ! gv(+1, 0, 0) = extrap3(f,i1p1,i2,j3,-1,0,0)
316  extrapWithMask(gv(+1, 0, 0), f, i1p1,i2,i3, i1p1,i2,j3, -1,0,0 )
317  else
318  gv(+1, 0, 0) = f(i1p1,i2,j3)
319  end if
320  ! grr=FRR(i1,i2,j3)
321  grr = GVRR()
322 
323  i2m1 = i2-1
324  if( i2m1.lt.mm2a .or. mask(i1,i2m1,i3).le.0 )then
325  ! f(i1,i2m1,j3)= extrap3(f,i1,i2m1,j3, 0,1,0)
326  ! gv( 0,-1, 0) = extrap3(f,i1,i2m1,j3, 0,1,0)
327  extrapWithMask(gv( 0,-1, 0), f, i1,i2m1,i3, i1,i2m1,j3, 0,1,0 )
328  else
329  gv( 0,-1, 0) = f(i1,i2m1,j3)
330  end if
331  i2p1 = i2+1
332  if( i2p1.gt.mm2b .or. mask(i1,i2p1,i3).le.0 )then
333  ! f(i1,i2p1,j3)= extrap3(f,i1,i2p1,j3, 0,-1,0)
334  ! gv( 0,+1, 0) = extrap3(f,i1,i2p1,j3, 0,-1,0)
335  extrapWithMask(gv( 0,+1, 0), f, i1,i2p1,i3, i1,i2p1,j3, 0,-1,0 )
336  else
337  gv( 0,+1, 0) = f(i1,i2p1,j3)
338  end if
339  ! gss=FSS(i1,i2,j3)
340  gss = GVSS()
341 
342  ! if( i1.eq.mm1a )then
343  ! grr =FRRa(i1,i2,j3)
344  ! else if( i1.eq.mm1b )then
345  ! grr =FRRb(i1,i2,j3)
346  ! else
347  ! grr=FRR(i1,i2,j3)
348  ! end if
349  ! if( i2.eq.mm2a )then
350  ! gss =FSSa(i1,i2,j3)
351  ! else if( i2.eq.mm2b )then
352  ! gss =FSSb(i1,i2,j3)
353  ! else
354  ! gss=FSS(i1,i2,j3)
355  ! end if
356 
357  #Else
358  stop 7
359  #End
360 
361  #Else
362  g=0.
363  ff=0.
364  #If #DIR eq "R"
365  ffr=0.
366  gs=0.
367  gss=0.
368  #Else
369  ffs=0.
370  gr=0.
371  grr=0.
372  #End
373 
374  #End
375 
376 #Elif #GRIDTYPE eq "curvilinear"
377 
378 ! Curvilinear grids use dr:
379 #defineMacro D12(axis) d12(axis)
380 #defineMacro D22(axis) d22(axis)
381 
382  #If #FORCING eq "forcing"
383  g = f(j1,j2,j3)
384  ff= f(i1,i2,i3)
385 
386  #If #DIM eq "3"
387  ax1 = mod(axis+1,nd)
388  ax2 = mod(axis+2,nd)
389  mdim(0,0)=mm1a
390  mdim(1,0)=mm1b
391  mdim(0,1)=mm2a
392  mdim(1,1)=mm2b
393  mdim(0,2)=mm3a
394  mdim(1,2)=mm3b
395  #End
396 
397 
398 
399  #If #DIR eq "R"
400  ! 2nd-order one sided:
401  ! ffr=is1*(-f(i1+2*is1,i2,i3)+4.*f(i1+is1,i2,i3)-3.*ff)*d12(0)
402  ! 3rd-order one sided:
403  ffr=is1*(-11.*ff+18.*f(i1+is1,i2,i3)-9.*f(i1+2*is1,i2,i3)+2.*f(i1+3*is1,i2,i3))/(6.*dr(0))
404 
405  ! NOTE: the forcing f and g are only assumed to be given where mask>0
406  ! In order to compute tangential derivatives of the forcing we may need to fill in
407  ! neighbouring values of the forcing at interp and unused points
408  fv( 0, 0, 0) = f(i1,i2,i3)
409  gv( 0, 0, 0) = f(j1,i2,i3)
410  i2m1 = i2-1
411  if( i2m1.lt.mm2a .or. mask(i1,i2m1,i3).le.0 )then
412  ! NOTE: We DO need to extrap f and g
413  ! f(i1,i2m1,i3)= extrap3(f,i1,i2m1,i3, 0,1,0)
414  ! f(j1,i2m1,i3)= extrap3(f,j1,i2m1,i3, 0,1,0)
415  ! fv( 0,-1, 0) = extrap3(f,i1,i2m1,i3, 0,1,0)
416  ! gv( 0,-1, 0) = extrap3(f,j1,i2m1,i3, 0,1,0)
417 
418  extrapWithMask( fv( 0,-1, 0), f, i1,i2m1,i3, i1,i2m1,i3, 0,1,0 )
419  extrapWithMask( gv( 0,-1, 0), f, i1,i2m1,i3, j1,i2m1,i3, 0,1,0 )
420 
421  else
422  fv( 0,-1, 0) = f(i1,i2m1,i3)
423  gv( 0,-1, 0) = f(j1,i2m1,i3)
424  end if
425  i2p1 = i2+1
426  if( i2p1.gt.mm2b .or. mask(i1,i2p1,i3).le.0 )then
427  ! f(i1,i2p1,i3)= extrap3(f,i1,i2p1,i3, 0,-1,0)
428  ! f(j1,i2p1,i3)= extrap3(f,j1,i2p1,i3, 0,-1,0)
429  ! fv( 0,+1, 0) = extrap3(f,i1,i2p1,i3, 0,-1,0)
430  ! gv( 0,+1, 0) = extrap3(f,j1,i2p1,i3, 0,-1,0)
431  extrapWithMask(fv( 0,+1, 0), f, i1,i2p1,i3, i1,i2p1,i3, 0,-1,0 )
432  extrapWithMask(gv( 0,+1, 0), f, i1,i2p1,i3, j1,i2p1,i3, 0,-1,0 )
433 
434  else
435  fv( 0,+1, 0) = f(i1,i2p1,i3)
436  gv( 0,+1, 0) = f(j1,i2p1,i3)
437  end if
438  ! ffs= FS(i1,i2,i3)
439  ! gs = FS(j1,i2,i3)
440  ! gss=FSS(j1,i2,i3)
441  ffs = FVS()
442  gs = GVS()
443  gss = GVSS()
444 
445  #If #DIM eq "3"
446 
447  i3m1 = i3-1
448  if( i3m1.lt.mm3a .or. mask(i1,i2,i3m1).le.0 )then
449  ! f(i1,i2,i3m1)= extrap3(f,i1,i2,i3m1, 0,0,1)
450  ! f(j1,i2,i3m1)= extrap3(f,j1,i2,i3m1, 0,0,1)
451  ! fv( 0, 0,-1) = extrap3(f,i1,i2,i3m1, 0,0,1)
452  ! gv( 0, 0,-1) = extrap3(f,j1,i2,i3m1, 0,0,1)
453  extrapWithMask(fv( 0, 0,-1), f, i1,i2,i3m1, i1,i2,i3m1, 0,0,1 )
454  extrapWithMask(gv( 0, 0,-1), f, i1,i2,i3m1, j1,i2,i3m1, 0,0,1 )
455 
456  else
457  fv( 0, 0,-1) = f(i1,i2,i3m1)
458  gv( 0, 0,-1) = f(j1,i2,i3m1)
459  end if
460  i3p1 = i3+1
461  if( i3p1.gt.mm3b .or. mask(i1,i2,i3p1).le.0 )then
462  ! f(i1,i2,i3p1)= extrap3(f,i1,i2,i3p1, 0,0,-1)
463  ! f(j1,i2,i3p1)= extrap3(f,j1,i2,i3p1, 0,0,-1)
464  ! fv( 0, 0,+1) = extrap3(f,i1,i2,i3p1, 0,0,-1)
465  ! gv( 0, 0,+1) = extrap3(f,j1,i2,i3p1, 0,0,-1)
466  extrapWithMask(fv( 0, 0,+1), f, i1,i2,i3p1, i1,i2,i3p1, 0,0,-1 )
467  extrapWithMask(gv( 0, 0,+1), f, i1,i2,i3p1, j1,i2,i3p1, 0,0,-1 )
468  else
469  fv( 0, 0,+1) = f(i1,i2,i3p1)
470  gv( 0, 0,+1) = f(j1,i2,i3p1)
471  end if
472  ! fft= FT(i1,i2,i3)
473  ! gt = FT(j1,i2,i3)
474  ! gtt=FTT(j1,i2,i3)
475  fft = FVT()
476  gt = GVT()
477  gtt = GVTT()
478 
479  ! compute the cross derivative: gst
480  ! Near physical or interpolation boundaries we may need to use a one sided approximation
481 
482  ! Evaluate g at neighbouring points so we can evaluate the cross derivative
483  getAdjacentBoundaryForcing(j1,i2m1,i3m1, j1,i2,i3, i1,i2m1,i3m1)
484  getAdjacentBoundaryForcing(j1,i2p1,i3m1, j1,i2,i3, i1,i2p1,i3m1)
485  getAdjacentBoundaryForcing(j1,i2m1,i3p1, j1,i2,i3, i1,i2m1,i3p1)
486  getAdjacentBoundaryForcing(j1,i2p1,i3p1, j1,i2,i3, i1,i2p1,i3p1)
487  gst = GVST()
488 
489  #End
490 
491  #Elif #DIR eq "S"
492  ! 2nd-order one sided:
493  ! is2*ffs=(-f(i1,i2+2*is2,i3)+4.*f(i1,i2+is2,i3)-3.*ff)*d12(1)
494  ! 3rd-order one sided:
495  ffs=is2*(-11.*ff+18.*f(i1,i2+is2,i3)-9.*f(i1,i2+2*is2,i3)+2.*f(i1,i2+3*is2,i3))/(6.*dr(1))
496 
497  fv( 0, 0, 0) = f(i1,i2,i3)
498  gv( 0, 0, 0) = f(i1,j2,i3)
499 
500  i1m1 = i1-1
501  if( i1m1.lt.mm1a .or. mask(i1m1,i2,i3).le.0 )then
502  ! f(i1m1,i2,i3)= extrap3(f,i1m1,i2,i3, 1,0,0)
503  ! f(i1m1,j2,i3)= extrap3(f,i1m1,j2,i3, 1,0,0)
504  ! fv(-1, 0, 0) = extrap3(f,i1m1,i2,i3, 1,0,0)
505  ! gv(-1, 0, 0) = extrap3(f,i1m1,j2,i3, 1,0,0)
506  extrapWithMask(fv(-1, 0, 0), f, i1m1,i2,i3, i1m1,i2,i3, 1,0,0 )
507  extrapWithMask(gv(-1, 0, 0), f, i1m1,i2,i3, i1m1,j2,i3, 1,0,0 )
508  else
509  fv(-1, 0, 0) = f(i1m1,i2,i3)
510  gv(-1, 0, 0) = f(i1m1,j2,i3)
511  end if
512  i1p1 = i1+1
513  if( i1p1.gt.mm1b .or. mask(i1p1,i2,i3).le.0 )then
514  ! f(i1p1,i2,i3)= extrap3(f,i1p1,i2,i3,-1,0,0)
515  ! f(i1p1,j2,i3)= extrap3(f,i1p1,j2,i3,-1,0,0)
516  ! fv(+1, 0, 0) = extrap3(f,i1p1,i2,i3,-1,0,0)
517  ! gv(+1, 0, 0) = extrap3(f,i1p1,j2,i3,-1,0,0)
518  extrapWithMask(fv(+1, 0, 0), f, i1p1,i2,i3, i1p1,i2,i3, -1,0,0 )
519  extrapWithMask(gv(+1, 0, 0), f, i1p1,i2,i3, i1p1,j2,i3, -1,0,0 )
520  else
521  fv(+1, 0, 0) = f(i1p1,i2,i3)
522  gv(+1, 0, 0) = f(i1p1,j2,i3)
523  end if
524  ! ffr= FR(i1,i2,i3)
525  ! gr = FR(i1,j2,i3)
526  ! grr=FRR(i1,j2,i3)
527  ffr = FVR()
528  gr = GVR()
529  grr = GVRR()
530 
531  #If #DIM eq "3"
532 
533  i3m1 = i3-1
534  if( i3m1.lt.mm3a .or. mask(i1,i2,i3m1).le.0 )then
535  ! f(i1,i2,i3m1)= extrap3(f,i1,i2,i3m1, 0,0, 1)
536  ! f(i1,j2,i3m1)= extrap3(f,i1,j2,i3m1, 0,0, 1)
537  ! fv( 0, 0,-1) = extrap3(f,i1,i2,i3m1, 0,0, 1)
538  ! gv( 0, 0,-1) = extrap3(f,i1,j2,i3m1, 0,0, 1)
539  extrapWithMask(fv( 0, 0,-1), f, i1,i2,i3m1, i1,i2,i3m1, 0,0,1 )
540  extrapWithMask(gv( 0, 0,-1), f, i1,i2,i3m1, i1,j2,i3m1, 0,0,1 )
541  else
542  fv( 0, 0,-1) = f(i1,i2,i3m1)
543  gv( 0, 0,-1) = f(i1,j2,i3m1)
544  end if
545  i3p1 = i3+1
546  if( i3p1.gt.mm3b .or. mask(i1,i2,i3p1).le.0 )then
547  ! f(i1,i2,i3p1)= extrap3(f,i1,i2,i3p1, 0,0,-1)
548  ! f(i1,j2,i3p1)= extrap3(f,i1,j2,i3p1, 0,0,-1)
549  ! fv( 0, 0,+1) = extrap3(f,i1,i2,i3p1, 0,0,-1)
550  ! gv( 0, 0,+1) = extrap3(f,i1,j2,i3p1, 0,0,-1)
551  extrapWithMask(fv( 0, 0,+1), f, i1,i2,i3p1, i1,i2,i3p1, 0,0,-1 )
552  extrapWithMask(gv( 0, 0,+1), f, i1,i2,i3p1, i1,j2,i3p1, 0,0,-1 )
553  else
554  fv( 0, 0,+1) = f(i1,i2,i3p1)
555  gv( 0, 0,+1) = f(i1,j2,i3p1)
556  end if
557  ! fft= FT(i1,i2,i3)
558  ! gt = FT(i1,j2,i3)
559  ! gtt=FTT(i1,j2,i3)
560  fft = FVT()
561  gt = GVT()
562  gtt = GVTT()
563 
564  ! Evaluate g at neighbouring points so we can evaluate the cross derivative
565  getAdjacentBoundaryForcing(i1m1,j2,i3m1, i1,j2,i3, i1m1,i2,i3m1)
566  getAdjacentBoundaryForcing(i1p1,j2,i3m1, i1,j2,i3, i1p1,i2,i3m1)
567  getAdjacentBoundaryForcing(i1m1,j2,i3p1, i1,j2,i3, i1m1,i2,i3p1)
568  getAdjacentBoundaryForcing(i1p1,j2,i3p1, i1,j2,i3, i1p1,i2,i3p1)
569 
570  grt = GVRT()
571 
572  #End
573 
574 
575  #Elif #DIR eq "T"
576 
577  ! 3rd-order one sided:
578  fft=is3*(-11.*ff+18.*f(i1,i2,i3+is3)-9.*f(i1,i2,i3+2*is3)+2.*f(i1,i2,i3+3*is3))/(6.*dr(2))
579 
580  fv( 0, 0, 0) = f(i1,i2,i3)
581  gv( 0, 0, 0) = f(i1,i2,j3)
582 
583  i1m1 = i1-1
584  if( i1m1.lt.mm1a .or. mask(i1m1,i2,i3).le.0 )then
585  ! f(i1m1,i2,i3)= extrap3(f,i1m1,i2,i3, 1,0,0)
586  ! f(i1m1,i2,j3)= extrap3(f,i1m1,i2,j3, 1,0,0)
587  ! fv(-1, 0, 0) = extrap3(f,i1m1,i2,i3, 1,0,0)
588  ! gv(-1, 0, 0) = extrap3(f,i1m1,i2,j3, 1,0,0)
589  extrapWithMask(fv(-1, 0, 0), f, i1m1,i2,i3, i1m1,i2,i3, 1,0,0 )
590  extrapWithMask(gv(-1, 0, 0), f, i1m1,i2,i3, i1m1,i2,j3, 1,0,0 )
591  else
592  fv(-1, 0, 0) = f(i1m1,i2,i3)
593  gv(-1, 0, 0) = f(i1m1,i2,j3)
594  endif
595  i1p1 = i1+1
596  if( i1p1.gt.mm1b .or. mask(i1p1,i2,i3).le.0 )then
597  ! f(i1p1,i2,i3)= extrap3(f,i1p1,i2,i3,-1,0,0)
598  ! f(i1p1,i2,j3)= extrap3(f,i1p1,i2,j3,-1,0,0)
599  ! fv(+1, 0, 0) = extrap3(f,i1p1,i2,i3,-1,0,0)
600  ! gv(+1, 0, 0) = extrap3(f,i1p1,i2,j3,-1,0,0)
601  extrapWithMask(fv(+1, 0, 0), f, i1p1,i2,i3, i1p1,i2,i3, -1,0,0 )
602  extrapWithMask(gv(+1, 0, 0), f, i1p1,i2,i3, i1p1,i2,j3, -1,0,0 )
603  else
604  fv(+1, 0, 0) = f(i1p1,i2,i3)
605  gv(+1, 0, 0) = f(i1p1,i2,j3)
606  endif
607  ! ffr= FR(i1,i2,i3)
608  ! gr = FR(i1,i2,j3)
609  ! grr=FRR(i1,i2,j3)
610  ffr = FVR()
611  gr = GVR()
612  grr = GVRR()
613 
614  i2m1 = i2-1
615  if( i2m1.lt.mm2a .or. mask(i1,i2m1,i3).le.0 )then
616  ! f(i1,i2m1,i3)= extrap3(f,i1,i2m1,i3, 0,1,0)
617  ! f(i1,i2m1,j3)= extrap3(f,i1,i2m1,j3, 0,1,0)
618  ! fv( 0,-1, 0) = extrap3(f,i1,i2m1,i3, 0,1,0)
619  ! gv( 0,-1, 0) = extrap3(f,i1,i2m1,j3, 0,1,0)
620  extrapWithMask(fv( 0,-1, 0), f, i1,i2m1,i3, i1,i2m1,i3, 0,1,0 )
621  extrapWithMask(gv( 0,-1, 0), f, i1,i2m1,i3, i1,i2m1,j3, 0,1,0 )
622  else
623  fv( 0,-1, 0) = f(i1,i2m1,i3)
624  gv( 0,-1, 0) = f(i1,i2m1,j3)
625  endif
626  i2p1 = i2+1
627  if( i2p1.gt.mm2b .or. mask(i1,i2p1,i3).le.0 )then
628  ! f(i1,i2p1,i3)= extrap3(f,i1,i2p1,i3, 0,-1,0)
629  ! f(i1,i2p1,j3)= extrap3(f,i1,i2p1,j3, 0,-1,0)
630  ! fv( 0,+1, 0) = extrap3(f,i1,i2p1,i3, 0,-1,0)
631  ! gv( 0,+1, 0) = extrap3(f,i1,i2p1,j3, 0,-1,0)
632  extrapWithMask(fv( 0,+1, 0), f, i1,i2p1,i3, i1,i2p1,i3, 0,-1,0 )
633  extrapWithMask(gv( 0,+1, 0), f, i1,i2p1,i3, i1,i2p1,j3, 0,-1,0 )
634  else
635  fv( 0,+1, 0) = f(i1,i2p1,i3)
636  gv( 0,+1, 0) = f(i1,i2p1,j3)
637  endif
638  ! ffs= FS(i1,i2,i3)
639  ! gs = FS(i1,i2,j3)
640  ! gss=FSS(i1,i2,j3)
641  ffs = FVS()
642  gs = GVS()
643  gss = GVSS()
644 
645  ! Evaluate g at neighbouring points so we can evaluate the cross derivative
646  getAdjacentBoundaryForcing(i1m1,i2m1,j3, i1,i2,j3, i1m1,i2m1,i3 )
647  getAdjacentBoundaryForcing(i1p1,i2m1,j3, i1,i2,j3, i1p1,i2m1,i3 )
648  getAdjacentBoundaryForcing(i1m1,i2p1,j3, i1,i2,j3, i1m1,i2p1,i3 )
649  getAdjacentBoundaryForcing(i1p1,i2p1,j3, i1,i2,j3, i1p1,i2p1,i3 )
650  grs = GVRS()
651 
652  #Else
653  stop 48480
654  #End
655 
656  #Else
657  ff=0.
658  ffr=0.
659  ffs=0.
660  fft=0.
661  g=0.
662  #If #DIR eq "R"
663  gs=0.
664  gss=0.
665  gt=0.
666  gtt=0.
667  gst=0.
668  #Elif #DIR eq "S"
669  gr=0.
670  grr=0.
671  gt=0.
672  gtt=0.
673  grt=0.
674  #Elif #DIR eq "T"
675  gr=0.
676  grr=0.
677  gs=0.
678  gss=0.
679  grs=0.
680  #Else
681  stop 48481
682  #End
683 
684  #End
685 
686 #Else
687  write(*,*) "Ogmg:NeuEqn:ERROR unknown gridType: GRIDTYPE"
688  stop 1863
689 #End
690 
691 #endMacro
692 
Generated on Fri Jan 4 2013 10:17:56 for Overture by   doxygen 1.8.3