MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_geometry.t
Go to the documentation of this file.
1 !> Module with geometry-related routines (e.g., divergence, curl)
3  implicit none
4  public
5 
6  integer :: coordinate=-1
7  integer, parameter :: cartesian = 0
8  integer, parameter :: cartesian_stretched= 1
9  integer, parameter :: cylindrical = 2
10  integer, parameter :: spherical = 3
11  integer, parameter :: cartesian_expansion= 4
12 
13 contains
14 
15  !> Set the coordinate system to be used
16  subroutine set_coordinate_system(geom)
18  character(len=*), intent(in) :: geom !< Name of the coordinate system
19 
20  ! Store the geometry name
21  geometry_name = geom
22 
23  select case (geom)
24  case ("Cartesian","Cartesian_1D","Cartesian_2D","Cartesian_3D")
25  ndir = ndim
27  case ("Cartesian_1D_expansion")
28  if (ndim /= 1) call mpistop("Geometry Cartesian_1D_expansion but ndim /= 1")
29  ndir = ndim
31  case ("Cartesian_1.5D")
32  if (ndim /= 1) call mpistop("Geometry Cartesian_1.5D but ndim /= 1")
34  ndir = 2
35  case ("Cartesian_1.75D")
36  if (ndim /= 1) call mpistop("Geometry Cartesian_1.75D but ndim /= 1")
38  ndir = 3
39  case ("Cartesian_2.5D")
40  if (ndim /= 2) call mpistop("Geometry Cartesian_2.5D but ndim /= 2")
42  ndir = 3
43  case ("cylindrical","cylindrical_2D","cylindrical_3D")
44  ndir = ndim
45  r_ = 1
46  z_ = 2
47  if(ndir==3) phi_ = 3
49  case ("cylindrical_2.5D")
50  if (ndim /= 2) call mpistop("Geometry cylindrical_2.5D but ndim /= 2")
51  ndir = 3
52  r_ = 1
53  z_ = 2
54  phi_ = 3
56  case ("polar","polar_2D","polar_3D")
57  ndir = ndim
58  r_ = 1
59  phi_ = 2
60  if(ndir==3) z_ = 3
62  case ("polar_1.5D")
63  if (ndim /= 1) call mpistop("Geometry polar_1.5D but ndim /= 1")
64  ndir = 2
65  r_ = 1
66  phi_ = 2
68  case ("polar_2.5D")
69  if (ndim /= 2) call mpistop("Geometry polar_2.5D but ndim /= 2")
70  ndir = 3
71  r_ = 1
72  phi_ = 2
73  z_ = 3
75  case ("spherical","spherical_2D","spherical_3D")
76  ndir = ndim
77  r_ = 1
78  if(ndir==3) phi_ = 3
79  z_ = -1
81  case ("spherical_2.5D")
82  if (ndim /= 2) &
83  call mpistop("Geometry spherical_2.5D requires ndim == 2")
84  ndir = 3
85  r_ = 1
86  phi_ = 3
87  z_ = -1
89  case default
90  call mpistop("Unknown geometry specified")
91  end select
92  end subroutine set_coordinate_system
93 
94  subroutine set_pole
96 
97  select case (coordinate)
98  case (spherical) {^ifthreed
99  ! For spherical grid, check whether phi-direction is periodic
100  if(periodb(ndim)) then
101  if(phi_/=3) call mpistop("phi_ should be 3 in 3D spherical coord!")
102  if(mod(ng3(1),2)/=0) &
103  call mpistop("Number of meshes in phi-direction should be even!")
104  if(abs(xprobmin2)<smalldouble) then
105  if(mype==0) write(unitterm,*) &
106  "Will apply pi-periodic conditions at northpole!"
107  poleb(1,2)=.true.
108  else
109  if(mype==0) write(unitterm,*) "There is no northpole!"
110  end if
111  if(abs(xprobmax2-dpi)<smalldouble) then
112  if(mype==0) write(unitterm,*) &
113  "Will apply pi-periodic conditions at southpole!"
114  poleb(2,2)=.true.
115  else
116  if(mype==0) write(unitterm,*) "There is no southpole!"
117  end if
118  end if}
119  case (cylindrical)
120  {
121  if (^d == phi_ .and. periodb(^d)) then
122  if(mod(ng^d(1),2)/=0) then
123  call mpistop("Number of meshes in phi-direction should be even!")
124  end if
125 
126  if(abs(xprobmin1)<smalldouble) then
127  if (mype==0) then
128  write(unitterm,*) "Will apply pi-periodic conditions at r=0"
129  end if
130  poleb(1,1)=.true.
131  else
132  if (mype==0) then
133  write(unitterm,*) "There is no cylindrical axis!"
134  end if
135  end if
136  end if\}
137  end select
138 
139  end subroutine set_pole
140 
141  !> Deallocate geometry-related variables
142  subroutine putgridgeo(igrid)
144  integer, intent(in) :: igrid
145 
146  deallocate(ps(igrid)%surfaceC,ps(igrid)%surface,ps(igrid)%dvolume,ps(igrid)%dx,&
147  psc(igrid)%dx,ps(igrid)%ds,psc(igrid)%dvolume)
148 
149  end subroutine putgridgeo
150 
151  !> calculate area of surfaces of cells
152  subroutine get_surface_area(s,ixG^L)
155 
156  type(state) :: s
157  integer, intent(in) :: ixG^L
158 
159  double precision :: x(ixg^s,ndim), xext(ixgmin1-1:ixgmax1,1), drs(ixg^s), dx2(ixg^s), dx3(ixg^s)
160  double precision :: exp_factor(ixg^s),del_exp_factor(ixg^s),exp_factor_primitive(ixg^s)
161 
162  select case (coordinate)
163 
164  case (cartesian_expansion)
165  drs(ixg^s)=s%dx(ixg^s,1)
166  x(ixg^s,1)=s%x(ixg^s,1)
167  {^ifoned
168  if(associated(usr_set_surface)) call usr_set_surface(ixg^l,x,drs,exp_factor,del_exp_factor,exp_factor_primitive)
169  s%surface(ixg^s,1)=exp_factor(ixg^s)
170  xext(0,1)=x(1,1)-half*drs(1)
171  xext(ixg^s,1)=x(ixg^s,1)+half*drs(ixg^s)
172  if(associated(usr_set_surface)) call usr_set_surface(ixgmin1-1,ixgmax1,xext,drs,exp_factor,del_exp_factor,exp_factor_primitive)
173  if (any(exp_factor .le. zero)) call mpistop("The area must always be strictly positive!")
174  s%surfaceC(ixgmin1-1:ixgmax1,1)=exp_factor(ixgmin1-1:ixgmax1)
175  }
176 
178  drs(ixg^s)=s%dx(ixg^s,1)
179  {^nooned
180  dx2(ixg^s)=s%dx(ixg^s,2)}
181  {^ifthreed
182  dx3(ixg^s)=s%dx(ixg^s,3)}
183 
184  {^ifoned
185  s%surfaceC(ixg^s,1)=1.d0
186  s%surface(ixg^s,1) =1.d0
187  }
188  {^iftwod
189  s%surfaceC(ixg^s,1)=dx2(ixg^s)
190  s%surfaceC(ixg^s,2)=drs(ixg^s)
191  s%surface(ixg^s,1) =dx2(ixg^s)
192  s%surface(ixg^s,2)=drs(ixg^s)
193  }
194  {^ifthreed
195  s%surfaceC(ixg^s,1)= dx2(ixg^s)*dx3(ixg^s)
196  s%surfaceC(ixg^s,2)= drs(ixg^s)*dx3(ixg^s)
197  s%surfaceC(ixg^s,3)= drs(ixg^s)*dx2(ixg^s)
198  s%surface(ixg^s,1)=s%surfaceC(ixg^s,1)
199  s%surface(ixg^s,2)=s%surfaceC(ixg^s,2)
200  s%surface(ixg^s,3)=s%surfaceC(ixg^s,3)
201  }
202  {s%surfaceC(0^d%ixG^s,^d)=s%surfaceC(1^d%ixG^s,^d);\}
203  case (spherical)
204  x(ixg^s,1)=s%x(ixg^s,1)
205  {^nooned
206  x(ixg^s,2)=s%x(ixg^s,2)
207  }
208  drs(ixg^s)=s%dx(ixg^s,1)
209  {^nooned
210  dx2(ixg^s)=s%dx(ixg^s,2)
211  }
212  {^ifthreed
213  dx3(ixg^s)=s%dx(ixg^s,3)
214  }
215 
216  s%surfaceC(ixg^s,1)=(x(ixg^s,1)+half*drs(ixg^s))**2 {^nooned &
217  *two*dsin(x(ixg^s,2))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
218 
219  {^nooned
220  s%surfaceC(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
221  *dsin(x(ixg^s,2)+half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
222 
223  {^ifthreed
224  s%surfaceC(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
225  }
226 
227  {^ifoned
228  s%surfaceC(0,1)=dabs(x(1,1)-half*drs(1))**2
229  }
230  {^iftwod
231  s%surfaceC(0,ixgmin2:ixgmax2,1)=(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
232  ixgmin2:ixgmax2))**2*two*dsin(x(1,ixgmin2:ixgmax2,2))*dsin(half*dx2(1,&
233  ixgmin2:ixgmax2))
234  s%surfaceC(ixgmin1:ixgmax1,0,2)=x(ixgmin1:ixgmax1,1,&
235  1)*drs(ixgmin1:ixgmax1,1)*dsin(x(ixgmin1:ixgmax1,1,&
236  2)-half*dx2(ixgmin1:ixgmax1,1))
237  }
238  {^ifthreed
239  s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=(x(1,ixgmin2:ixgmax2,&
240  ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,&
241  ixgmin3:ixgmax3))**2*two*dsin(x(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3,&
242  2))*dsin(half*dx2(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx3(1,&
243  ixgmin2:ixgmax2,ixgmin3:ixgmax3)
244  s%surfaceC(ixgmin1:ixgmax1,0,ixgmin3:ixgmax3,2)=x(ixgmin1:ixgmax1,1,&
245  ixgmin3:ixgmax3,1)*drs(ixgmin1:ixgmax1,1,&
246  ixgmin3:ixgmax3)*dsin(x(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3,&
247  2)-half*dx2(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3))*dx3(ixgmin1:ixgmax1,1,&
248  ixgmin3:ixgmax3)
249  s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,0,3)=&
250  s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,1,3)
251  }
252 
253  s%surface(ixg^s,1)=x(ixg^s,1)**2 {^nooned &
254  *two*dsin(x(ixg^s,2))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
255  {^nooned
256  s%surface(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
257  *dsin(x(ixg^s,2))}{^ifthreed*dx3(ixg^s)}
258 
259  {^ifthreed
260  s%surface(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)}
261 
262  case (cylindrical)
263  x(ixg^s,1)=s%x(ixg^s,1)
264  drs(ixg^s)=s%dx(ixg^s,1)
265  {^nooned
266  dx2(ixg^s)=s%dx(ixg^s,2)}
267  {^ifthreed
268  dx3(ixg^s)=s%dx(ixg^s,3)}
269 
270  s%surfaceC(ixg^s,1)=dabs(x(ixg^s,1)+half*drs(ixg^s)){^de&*dx^de(ixg^s) }
271  {^nooned
272  if (z_==2) s%surfaceC(ixg^s,2)=x(ixg^s,1)*drs(ixg^s){^ifthreed*dx3(ixg^s)}
273  if (phi_==2) s%surfaceC(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}
274  }
275  {^ifthreed
276  if (z_==3) s%surfaceC(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
277  if (phi_==3) s%surfaceC(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)
278  }
279  {^ifoned
280  s%surfaceC(0,1)=dabs(x(1,1)-half*drs(1))
281  }
282  {^iftwod
283  s%surfaceC(0,ixgmin2:ixgmax2,1)=dabs(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
284  ixgmin2:ixgmax2))*dx2(1,ixgmin2:ixgmax2)
285  }
286  {^ifthreed
287  s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=dabs(x(1,ixgmin2:ixgmax2,&
288  ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx2(1,&
289  ixgmin2:ixgmax2,ixgmin3:ixgmax3)*dx3(1,ixgmin2:ixgmax2,&
290  ixgmin3:ixgmax3)
291  }
292  {s%surfaceC(0^de%ixG^s,^de)=s%surfaceC(1^de%ixG^s,^de);\}
293 
294  s%surface(ixg^s,1)=dabs(x(ixg^s,1)){^de&*dx^de(ixg^s) }
295  {^nooned
296  if (z_==2) s%surface(ixg^s,2)=x(ixg^s,1)*drs(ixg^s){^ifthreed*dx3(ixg^s)}
297  if (phi_==2) s%surface(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}}
298  {^ifthreed
299  if (z_==3) s%surface(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
300  if (phi_==3) s%surface(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)}
301 
302  case default
303  call mpistop("Sorry, coordinate unknown")
304  end select
305 
306  end subroutine get_surface_area
307 
308  !> Calculate gradient of a scalar q within ixL in direction idir
309  subroutine gradient(q,ixI^L,ixO^L,idir,gradq)
311 
312  integer, intent(in) :: ixI^L, ixO^L, idir
313  double precision, intent(in) :: q(ixi^s)
314  double precision, intent(inout) :: gradq(ixi^s)
315  double precision :: x(ixi^s,1:ndim)
316  integer :: jxO^L, hxO^L
317 
318  x(ixi^s,1:ndim)=block%x(ixi^s,1:ndim)
319 
320  hxo^l=ixo^l-kr(idir,^d);
321  jxo^l=ixo^l+kr(idir,^d);
322  select case(coordinate)
323  case(cartesian)
324  gradq(ixo^s)=half*(q(jxo^s)-q(hxo^s))/dxlevel(idir)
326  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
327  case(spherical)
328  select case(idir)
329  case(1)
330  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,1)-x(hxo^s,1)))
331  {^nooned
332  case(2)
333  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
334  }
335  {^ifthreed
336  case(3)
337  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,3)-x(hxo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)))
338  }
339  end select
340  case(cylindrical)
341  if(idir==phi_) then
342  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,phi_)-x(hxo^s,phi_))*x(ixo^s,r_))
343  else
344  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
345  end if
346  case default
347  call mpistop('Unknown geometry')
348  end select
349 
350  end subroutine gradient
351 
352  !> Calculate gradient of a scalar q in direction idir at cell interfaces
353  subroutine gradientx(q,x,ixI^L,ixO^L,idir,gradq,fourth_order)
355 
356  integer, intent(in) :: ixI^L, ixO^L, idir
357  double precision, intent(in) :: q(ixi^s), x(ixi^s,1:ndim)
358  double precision, intent(inout) :: gradq(ixi^s)
359  logical, intent(in) :: fourth_order
360  integer :: jxO^L, hxO^L, kxO^L
361 
362  if(fourth_order) then
363  ! Fourth order, stencil width is two
364  kxo^l=ixo^l^ladd2;
365  if(iximin^d>kxomin^d.or.iximax^d<kxomax^d|.or.) &
366  call mpistop("Error in gradientx: Non-conforming input limits")
367  hxo^l=ixo^l-kr(idir,^d);
368  jxo^l=ixo^l+kr(idir,^d);
369  kxo^l=ixo^l+kr(idir,^d)*2;
370  else
371  hxo^l=ixo^l;
372  end if
373  jxo^l=ixo^l+kr(idir,^d);
374  select case(coordinate)
375  case(cartesian)
376  if(fourth_order) then
377  gradq(ixo^s)=(27.d0*(q(jxo^s)-q(ixo^s))-q(kxo^s)+q(hxo^s))/24.d0/dxlevel(idir)
378  else
379  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/dxlevel(idir)
380  end if
382  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
383  case(spherical)
384  select case(idir)
385  case(1)
386  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,1)-x(hxo^s,1)))
387  {^nooned
388  case(2)
389  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
390  }
391  {^ifthreed
392  case(3)
393  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,3)-x(hxo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)))
394  }
395  end select
396  case(cylindrical)
397  if(idir==phi_) then
398  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,phi_)-x(hxo^s,phi_))*x(ixo^s,r_))
399  else
400  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
401  end if
402  case default
403  call mpistop('Unknown geometry')
404  end select
405 
406  end subroutine gradientx
407 
408  !> Calculate gradient of a scalar q within ixL in direction idir
409  !> first use limiter to go from cell center to edge
410  subroutine gradients(q,ixI^L,ixO^L,idir,gradq)
412  use mod_limiter
413  use mod_ppm
414 
415  integer, intent(in) :: ixI^L, ixO^L, idir
416  double precision, intent(in) :: q(ixi^s)
417  double precision, intent(inout) :: gradq(ixi^s)
418  double precision ,dimension(ixI^S) :: qC,qL,qR,dqC,ldq,rdq
419 
420  double precision :: x(ixi^s,1:ndim)
421  double precision :: invdx
422  integer :: hxO^L,ixC^L,jxC^L,gxC^L,hxC^L
423 
424  x(ixi^s,1:ndim)=block%x(ixi^s,1:ndim)
425 
426  invdx=1.d0/dxlevel(idir)
427  hxo^l=ixo^l-kr(idir,^d);
428  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
429  jxc^l=ixc^l+kr(idir,^d);
430  gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
431  hxc^l=gxc^l+kr(idir,^d);
432 
433  ! set the gradient limiter here
434  qr(gxc^s) = q(hxc^s)
435  ql(gxc^s) = q(gxc^s)
436  if (typegradlimiter/=limiter_ppm) then
437  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
438  call dwlimiter2(dqc,ixi^l,gxc^l,idir,typegradlimiter,ldw=ldq,rdw=rdq)
439  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
440  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
441  else
442  call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
443  endif
444 
445  select case(coordinate)
446  case(cartesian)
447  gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))*invdx
449  gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
450  case(spherical)
451  gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
452  select case(idir)
453  case(2)
454  gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
455  {^ifthreed
456  case(3)
457  gradq(ixo^s)=gradq(ixo^s)/(x(ixo^s,1)*dsin(x(ixo^s,2)))
458  }
459  end select
460  case(cylindrical)
461  gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
462  if(idir==phi_) gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
463  end select
464 
465  end subroutine gradients
466 
467  !> Calculate divergence of a vector qvec within ixL
468  subroutine divvector(qvec,ixI^L,ixO^L,divq,fourthorder,sixthorder)
470  integer, intent(in) :: ixI^L,ixO^L
471  double precision, intent(in) :: qvec(ixi^s,1:ndir)
472  double precision, intent(inout) :: divq(ixi^s)
473  logical, intent(in), optional :: fourthorder !< Default: false
474  logical, intent(in), optional :: sixthorder !< Default: false
475  logical :: use_4th_order
476  logical :: use_6th_order
477  double precision :: qC(ixi^s), invdx(1:ndim)
478  integer :: jxO^L, hxO^L, ixC^L, jxC^L
479  integer :: idims, ix^L, gxO^L, kxO^L
480  integer :: lxO^L, fxO^L
481 
482  use_4th_order = .false.
483  use_6th_order = .false.
484  if (present(fourthorder)) use_4th_order = fourthorder
485  if (present(sixthorder)) use_6th_order = sixthorder
486  if(use_4th_order .and. use_6th_order) &
487  call mpistop("divvector: using 4th and 6th order at the same time")
488 
489  if(use_4th_order) then
490  if (.not. slab_uniform) &
491  call mpistop("divvector: 4th order only supported for slab geometry")
492  ! Fourth order, stencil width is two
493  ix^l=ixo^l^ladd2;
494  else if(use_6th_order) then
495  ! Sixth order, stencil width is three
496  if (.not. slab_uniform) &
497  call mpistop("divvector: 6th order only supported for slab geometry")
498  ix^l=ixo^l^ladd3;
499  else
500  ! Second order, stencil width is one
501  ix^l=ixo^l^ladd1;
502  end if
503 
504  if (iximin^d>ixmin^d.or.iximax^d<ixmax^d|.or.) &
505  call mpistop("Error in divvector: Non-conforming input limits")
506 
507  invdx=1.d0/dxlevel
508  divq(ixo^s)=0.0d0
509 
510  if (slab_uniform) then
511  do idims=1,ndim
512  if(use_6th_order) then
513  lxo^l=ixo^l+3*kr(idims,^d);
514  kxo^l=ixo^l+2*kr(idims,^d);
515  jxo^l=ixo^l+kr(idims,^d);
516  hxo^l=ixo^l-kr(idims,^d);
517  gxo^l=ixo^l-2*kr(idims,^d);
518  fxo^l=ixo^l-3*kr(idims,^d);
519  divq(ixo^s)=divq(ixo^s)+&
520  (-qvec(fxo^s,idims)+9.d0*qvec(gxo^s,idims)-45.d0*qvec(hxo^s,idims)&
521  +qvec(lxo^s,idims)-9.d0*qvec(kxo^s,idims)+45.d0*qvec(jxo^s,idims))&
522  /(60.d0*dxlevel(idims))
523  else if(use_4th_order) then
524  ! Use fourth order scheme
525  kxo^l=ixo^l+2*kr(idims,^d);
526  jxo^l=ixo^l+kr(idims,^d);
527  hxo^l=ixo^l-kr(idims,^d);
528  gxo^l=ixo^l-2*kr(idims,^d);
529  divq(ixo^s)=divq(ixo^s)+&
530  (-qvec(kxo^s,idims)+8.d0*qvec(jxo^s,idims)-8.d0*&
531  qvec(hxo^s,idims)+qvec(gxo^s,idims))/(12.d0*dxlevel(idims))
532  else
533  ! Use second order scheme
534  jxo^l=ixo^l+kr(idims,^d);
535  hxo^l=ixo^l-kr(idims,^d);
536  divq(ixo^s)=divq(ixo^s)+half*(qvec(jxo^s,idims) &
537  - qvec(hxo^s,idims))*invdx(idims)
538  end if
539  end do
540  else
541  do idims=1,ndim
542  hxo^l=ixo^l-kr(idims,^d);
543  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
544  jxc^l=ixc^l+kr(idims,^d);
545  if(stretched_dim(idims) .and. stretch_uncentered) then
546  ! linear interpolation at cell interface along stretched dimension
547  qc(ixc^s)=block%surfaceC(ixc^s,idims)*(qvec(ixc^s,idims)+0.5d0*block%dx(ixc^s,idims)*&
548  (qvec(jxc^s,idims)-qvec(ixc^s,idims))/(block%x(jxc^s,idims)-block%x(ixc^s,idims)))
549  else
550  qc(ixc^s)=block%surfaceC(ixc^s,idims)*half*(qvec(ixc^s,idims)+qvec(jxc^s,idims))
551  end if
552  divq(ixo^s)=divq(ixo^s)+qc(ixo^s)-qc(hxo^s)
553  end do
554  divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
555  end if
556  end subroutine divvector
557 
558  !> Calculate divergence of a vector qvec within ixL
559  !> using limited extrapolation to cell edges
560  subroutine divvectors(qvec,ixI^L,ixO^L,divq)
562  use mod_limiter
563  use mod_ppm
564 
565  integer, intent(in) :: ixI^L,ixO^L
566  double precision, intent(in) :: qvec(ixi^s,1:ndir)
567  double precision, intent(inout) :: divq(ixi^s)
568  double precision, dimension(ixI^S) :: qL,qR,dqC,ldq,rdq
569 
570  double precision :: invdx(1:ndim)
571  integer :: hxO^L,ixC^L,jxC^L,idims,ix^L,gxC^L,hxC^L
572 
573  ix^l=ixo^l^ladd2;
574 
575  if (iximin^d>ixmin^d.or.iximax^d<ixmax^d|.or.) &
576  call mpistop("Error in divvectorS: Non-conforming input limits")
577 
578  invdx=1.d0/dxlevel
579  divq(ixo^s)=zero
580  do idims=1,ndim
581  hxo^l=ixo^l-kr(idims,^d);
582  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
583  jxc^l=ixc^l+kr(idims,^d);
584  gxcmin^d=ixcmin^d-kr(idims,^d);gxcmax^d=jxcmax^d;
585  hxc^l=gxc^l+kr(idims,^d);
586  qr(gxc^s) = qvec(hxc^s,idims)
587  ql(gxc^s) = qvec(gxc^s,idims)
588  if(typegradlimiter/=limiter_ppm) then
589  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
590  call dwlimiter2(dqc,ixi^l,gxc^l,idims,typegradlimiter,ldw=ldq,rdw=rdq)
591  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
592  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
593  else
594  dqc(ixi^s)=qvec(ixi^s,idims)
595  call ppmlimitervar(ixi^l,ixo^l,idims,dqc,dqc,ql,qr)
596  endif
597 
598  if (slab_uniform) then
599  divq(ixo^s)=divq(ixo^s)+half*(qr(ixo^s)-ql(hxo^s))*invdx(idims)
600  else
601  qr(ixc^s)=block%surfaceC(ixc^s,idims)*qr(ixc^s)
602  ql(ixc^s)=block%surfaceC(ixc^s,idims)*ql(ixc^s)
603  divq(ixo^s)=divq(ixo^s)+qr(ixo^s)-ql(hxo^s)
604  end if
605  end do
606  if(.not.slab_uniform) divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
607 
608  end subroutine divvectors
609 
610  !> Calculate curl of a vector qvec within ixL
611  !> Options to
612  !> employ standard second order CD evaluations
613  !> use Gauss theorem for non-Cartesian grids
614  !> use Stokes theorem for non-Cartesian grids
615  subroutine curlvector(qvec,ixI^L,ixO^L,curlvec,idirmin,idirmin0,ndir0,fourthorder)
617 
618  integer, intent(in) :: ixI^L,ixO^L
619  integer, intent(in) :: ndir0, idirmin0
620  integer, intent(inout) :: idirmin
621  double precision, intent(in) :: qvec(ixi^s,1:ndir0)
622  double precision, intent(inout) :: curlvec(ixi^s,idirmin0:3)
623  logical, intent(in), optional :: fourthorder !< Default: false
624 
625  integer :: ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L,kxO^L,gxO^L
626  double precision :: invdx(1:ndim)
627  double precision :: tmp(ixi^s),tmp2(ixi^s),xC(ixi^s),surface(ixi^s)
628  logical :: use_4th_order
629 
630  ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k
631  ! Curl can have components (idirmin:3)
632  ! Determine exact value of idirmin while doing the loop.
633 
634  use_4th_order = .false.
635  if (present(fourthorder)) use_4th_order = fourthorder
636 
637  if (use_4th_order) then
638  if (.not. slab_uniform) &
639  call mpistop("divvector: 4th order only supported for slab geometry")
640  ! Fourth order, stencil width is two
641  ixc^l=ixo^l^ladd2;
642  else
643  ! Second order, stencil width is one
644  ixc^l=ixo^l^ladd1;
645  end if
646 
647  if (iximin^d>ixcmin^d.or.iximax^d<ixcmax^d|.or.) &
648  call mpistop("Error in divvector: Non-conforming input limits")
649 
650  idirmin=4
651  curlvec(ixo^s,idirmin0:3)=zero
652 
653  ! all non-Cartesian cases
654  select case(coordinate)
655  case(cartesian) ! Cartesian grids
656  invdx=1.d0/dxlevel
657  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
658  if(lvc(idir,jdir,kdir)/=0)then
659  if (.not. use_4th_order) then
660  ! Use second order scheme
661  jxo^l=ixo^l+kr(jdir,^d);
662  hxo^l=ixo^l-kr(jdir,^d);
663  tmp(ixo^s)=half*(qvec(jxo^s,kdir) &
664  - qvec(hxo^s,kdir))*invdx(jdir)
665  else
666  ! Use fourth order scheme
667  kxo^l=ixo^l+2*kr(jdir,^d);
668  jxo^l=ixo^l+kr(jdir,^d);
669  hxo^l=ixo^l-kr(jdir,^d);
670  gxo^l=ixo^l-2*kr(jdir,^d);
671  tmp(ixo^s)=(-qvec(kxo^s,kdir) + 8.0d0 * qvec(jxo^s,kdir) - 8.0d0 * &
672  qvec(hxo^s,kdir) + qvec(gxo^s,kdir))/(12.0d0 * dxlevel(jdir))
673  end if
674  if(lvc(idir,jdir,kdir)==1)then
675  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
676  else
677  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
678  endif
679  if(idir<idirmin)idirmin=idir
680  endif
681  enddo; enddo; enddo;
682  case(cartesian_stretched) ! stretched Cartesian grids
683  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
684  if(lvc(idir,jdir,kdir)/=0)then
685  select case(typecurl)
686  case('central')
687  tmp(ixi^s)=qvec(ixi^s,kdir)
688  hxo^l=ixo^l-kr(jdir,^d);
689  jxo^l=ixo^l+kr(jdir,^d);
690  ! second order centered differencing
691  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,jdir)-block%x(hxo^s,jdir))
692  case('Gaussbased')
693  hxo^l=ixo^l-kr(jdir,^d);
694  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
695  jxc^l=ixc^l+kr(jdir,^d);
696  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
697  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
698  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
699  case('Stokesbased')
700  hxo^l=ixo^l-kr(jdir,^d);
701  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
702  jxc^l=ixc^l+kr(jdir,^d);
703  if(kdir<=ndim)then
704  tmp(ixc^s)=block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
705  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
706  else
707  tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
708  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
709  endif
710  if(idir<=ndim)then
711  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%surface(ixo^s,idir)
712  else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
713  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(block%ds(ixo^s,jdir)*block%ds(ixo^s,kdir))
714  endif
715  case default
716  call mpistop('no such curl evaluator')
717  end select
718  if(lvc(idir,jdir,kdir)==1)then
719  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
720  else
721  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
722  endif
723  if(idir<idirmin)idirmin=idir
724  endif
725  enddo; enddo; enddo;
726  case(spherical) ! possibly stretched spherical grids
727  select case(typecurl)
728  case('central') ! ok for any dimensionality
729  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
730  if(lvc(idir,jdir,kdir)/=0)then
731  tmp(ixi^s)=qvec(ixi^s,kdir)
732  hxo^l=ixo^l-kr(jdir,^d);
733  jxo^l=ixo^l+kr(jdir,^d);
734  select case(jdir)
735  case(1)
736  tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1)
737  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1))*block%x(ixo^s,1))
738  {^nooned case(2)
739  if(idir==1) tmp(ixi^s)=tmp(ixi^s)*dsin(block%x(ixi^s,2))
740  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
741  if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(block%x(ixo^s,2))
742  }
743  {^ifthreed case(3)
744  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1)*dsin(block%x(ixo^s,2)))
745  }
746  end select
747  if(lvc(idir,jdir,kdir)==1)then
748  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
749  else
750  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
751  endif
752  if(idir<idirmin)idirmin=idir
753  endif
754  enddo; enddo; enddo;
755  case('Gaussbased')
756  do idir=idirmin0,3;
757  do jdir=1,ndim; do kdir=1,ndir0
758  if(lvc(idir,jdir,kdir)/=0)then
759  hxo^l=ixo^l-kr(jdir,^d);
760  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
761  jxc^l=ixc^l+kr(jdir,^d);
762  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
763  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
764  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
765  if(lvc(idir,jdir,kdir)==1)then
766  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
767  else
768  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
769  endif
770  if(idir<idirmin)idirmin=idir
771  endif
772  enddo; enddo;
773  ! geometric terms
774  if(idir==2.and.phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,phi_)/block%x(ixo^s,r_)
775  {^nooned
776  if(idir==phi_) curlvec(ixo^s,phi_)=curlvec(ixo^s,phi_)-qvec(ixo^s,2)/block%x(ixo^s,r_) &
777  +qvec(ixo^s,r_)*dcos(block%x(ixo^s,2))/(block%x(ixo^s,r_)*dsin(block%x(ixo^s,2)))
778  }
779  enddo;
780  case('Stokesbased')
781  !if(ndim<3) call mpistop("Stokesbased for 3D spherical only")
782  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
783  if(lvc(idir,jdir,kdir)/=0)then
784  select case(idir)
785  case(1)
786  if(jdir<kdir) then
787  ! idir=1,jdir=2,kdir=3
788  !! integral along 3rd dimension
789  hxo^l=ixo^l-kr(jdir,^d);
790  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
791  jxc^l=ixc^l+kr(jdir,^d);
792  ! qvec(3) at cell interface along 2nd dimension
793  if(stretched_dim(jdir) .and. stretch_uncentered) then
794  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
795  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
796  else
797  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
798  end if
799  ! 2nd coordinate at cell interface along 2nd dimension
800  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
801  curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
802  dsin(xc(hxo^s))*tmp(hxo^s))*block%dx(ixo^s,kdir)
803  !! integral along 2nd dimension
804  hxo^l=ixo^l-kr(kdir,^d);
805  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
806  jxc^l=ixc^l+kr(kdir,^d);
807  ! qvec(2) at cell interface along 3rd dimension
808  if(stretched_dim(kdir) .and. stretch_uncentered) then
809  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
810  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
811  else
812  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
813  end if
814  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir))&
815  /block%surface(ixo^s,idir)*block%x(ixo^s,idir)
816  end if
817  case(2)
818  if(jdir<kdir) then
819  ! idir=2,jdir=1,kdir=3
820  !! integral along 1st dimension
821  hxo^l=ixo^l-kr(kdir,^d);
822  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
823  jxc^l=ixc^l+kr(kdir,^d);
824  ! qvec(1) at cell interface along 3rd dimension
825  if(stretched_dim(kdir) .and. stretch_uncentered) then
826  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
827  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
828  else
829  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
830  end if
831  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,1)
832  !! integral along 3rd dimension
833  hxo^l=ixo^l-kr(jdir,^d);
834  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
835  jxc^l=ixc^l+kr(jdir,^d);
836  ! qvec(3) at cell interface along 1st dimension
837  if(stretched_dim(jdir) .and. stretch_uncentered) then
838  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
839  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
840  else
841  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
842  end if
843  ! 1st coordinate at cell interface along 1st dimension
844  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
845  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
846  dsin(block%x(ixo^s,idir))*block%dx(ixo^s,kdir))/block%surface(ixo^s,idir)
847  end if
848  case(3)
849  if(jdir<kdir) then
850  ! idir=3,jdir=1,kdir=2
851  !! integral along 1st dimension
852  hxo^l=ixo^l-kr(kdir,^d);
853  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
854  jxc^l=ixc^l+kr(kdir,^d);
855  ! qvec(1) at cell interface along 2nd dimension
856  if(stretched_dim(kdir) .and. stretch_uncentered) then
857  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
858  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
859  else
860  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
861  end if
862  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
863  !! integral along 2nd dimension
864  hxo^l=ixo^l-kr(jdir,^d);
865  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
866  jxc^l=ixc^l+kr(jdir,^d);
867  ! qvec(2) at cell interface along 1st dimension
868  if(stretched_dim(jdir) .and. stretch_uncentered) then
869  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
870  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
871  else
872  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
873  end if
874  ! 1st coordinate at cell interface along 1st dimension
875  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
876  if(ndim==3) then
877  surface(ixo^s)=block%surface(ixo^s,idir)
878  else
879  surface(ixo^s)=block%x(ixo^s,jdir)*block%dx(ixo^s,kdir)*block%dx(ixo^s,jdir)
880  end if
881  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
882  /surface(ixo^s)
883  end if
884  end select
885  if(idir<idirmin)idirmin=idir
886  endif
887  enddo; enddo; enddo;
888  case default
889  call mpistop('no such curl evaluator')
890  end select
891  case(cylindrical) ! possibly stretched cylindrical grids
892  select case(typecurl)
893  case('central') ! works for any dimensionality, polar/cylindrical
894  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
895  if(lvc(idir,jdir,kdir)/=0)then
896  tmp(ixi^s)=qvec(ixi^s,kdir)
897  hxo^l=ixo^l-kr(jdir,^d);
898  jxo^l=ixo^l+kr(jdir,^d);
899  if(z_==3.or.z_==-1) then
900  ! Case Polar_2D, Polar_2.5D or Polar_3D, i.e. R,phi,Z
901  select case(jdir)
902  case(1)
903  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
904  ! computes d(R V_phi)/dR or d V_Z/dR
905  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
906  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
907  {^nooned case(2)
908  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
909  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
910  }
911  {^ifthreed case(3)
912  ! handles d V_phi/dZ or d V_R/dZ
913  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,3)-block%x(hxo^s,3))
914  }
915  end select
916  end if
917  if(phi_==3.or.phi_==-1) then
918  ! Case Cylindrical_2D, Cylindrical_2.5D or Cylindrical_3D, i.e. R,Z,phi
919  select case(jdir)
920  case(1)
921  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
922  ! computes d(R V_phi)/dR or d V_Z/dR
923  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
924  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
925  {^nooned case(2)
926  ! handles d V_phi/dZ or d V_R/dZ
927  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,2)-block%x(hxo^s,2))
928  }
929  {^ifthreed case(3)
930  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
931  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1))
932  }
933  end select
934  end if
935  if(lvc(idir,jdir,kdir)==1)then
936  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
937  else
938  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
939  endif
940  if(idir<idirmin)idirmin=idir
941  endif
942  enddo; enddo; enddo;
943  case('Gaussbased') ! works for any dimensionality, polar/cylindrical
944  if(ndim<2) call mpistop("Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
945  do idir=idirmin0,3;
946  do jdir=1,ndim; do kdir=1,ndir0
947  if(lvc(idir,jdir,kdir)/=0)then
948  hxo^l=ixo^l-kr(jdir,^d);
949  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
950  jxc^l=ixc^l+kr(jdir,^d);
951  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
952  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
953  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
954  if(lvc(idir,jdir,kdir)==1)then
955  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
956  else
957  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
958  endif
959  if(idir<idirmin)idirmin=idir
960  endif
961  enddo; enddo;
962  ! geometric term from d e_R/d phi= e_phi for unit vectors e_R, e_phi
963  ! but minus sign appears due to R,Z,phi ordering (?)
964  ! note that in cylindrical 2D (RZ), phi_ is -1
965  ! note that in polar 2D (Rphi), z_ is -1
966  if((idir==phi_.or.(phi_==-1.and.idir==3)).and.z_>0) then
967  ! cylindrical
968  if( z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,z_)/block%x(ixo^s,r_)
969  ! polar
970  if(phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,z_)/block%x(ixo^s,r_)
971  endif
972  enddo;
973  case('Stokesbased')
974  if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only")
975  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
976  if(lvc(idir,jdir,kdir)/=0)then
977  if(idir==r_) then
978  if(jdir==phi_) then
979  ! idir=r,jdir=phi,kdir=z
980  !! integral along z dimension
981  hxo^l=ixo^l-kr(jdir,^d);
982  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
983  jxc^l=ixc^l+kr(jdir,^d);
984  ! qvec(z) at cell interface along phi dimension
985  if(stretched_dim(jdir) .and. stretch_uncentered) then
986  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
987  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
988  else
989  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
990  end if
991  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,kdir)
992  !! integral along phi dimension
993  hxo^l=ixo^l-kr(kdir,^d);
994  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
995  jxc^l=ixc^l+kr(kdir,^d);
996  ! qvec(phi) at cell interface along z dimension
997  if(stretched_dim(kdir) .and. stretch_uncentered) then
998  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
999  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1000  else
1001  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1002  end if
1003  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%x(ixo^s,idir)*block%dx(ixo^s,jdir))&
1004  /block%surface(ixo^s,idir)
1005  end if
1006  else if(idir==phi_) then
1007  if(jdir<kdir) then
1008  ! idir=phi,jdir=r,kdir=z
1009  !! integral along r dimension
1010  hxo^l=ixo^l-kr(kdir,^d);
1011  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1012  jxc^l=ixc^l+kr(kdir,^d);
1013  ! qvec(r) at cell interface along z dimension
1014  if(stretched_dim(kdir) .and. stretch_uncentered) then
1015  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1016  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1017  else
1018  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1019  end if
1020  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,jdir)
1021  !! integral along z dimension
1022  hxo^l=ixo^l-kr(jdir,^d);
1023  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1024  jxc^l=ixc^l+kr(jdir,^d);
1025  ! qvec(z) at cell interface along r dimension
1026  if(stretched_dim(jdir) .and. stretch_uncentered) then
1027  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1028  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1029  else
1030  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1031  end if
1032  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,kdir))&
1033  /block%surface(ixo^s,idir)
1034  end if
1035  else ! idir==z_
1036  if(jdir<kdir) then
1037  ! idir=z,jdir=r,kdir=phi
1038  !! integral along r dimension
1039  hxo^l=ixo^l-kr(kdir,^d);
1040  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1041  jxc^l=ixc^l+kr(kdir,^d);
1042  ! qvec(r) at cell interface along phi dimension
1043  if(stretched_dim(kdir) .and. stretch_uncentered) then
1044  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1045  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1046  else
1047  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1048  end if
1049  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
1050  !! integral along phi dimension
1051  hxo^l=ixo^l-kr(jdir,^d);
1052  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1053  jxc^l=ixc^l+kr(jdir,^d);
1054  ! qvec(phi) at cell interface along r dimension
1055  if(stretched_dim(jdir) .and. stretch_uncentered) then
1056  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1057  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1058  else
1059  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1060  end if
1061  ! r coordinate at cell interface along r dimension
1062  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
1063  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
1064  /block%surface(ixo^s,idir)
1065  end if
1066  end if
1067  if(idir<idirmin)idirmin=idir
1068  endif
1069  enddo; enddo; enddo;
1070  case default
1071  call mpistop('no such curl evaluator')
1072  end select
1073  case default
1074  call mpistop('not possible to calculate curl')
1075  end select
1076 
1077  end subroutine curlvector
1078 
1079  !> Calculate idim transverse components of curl of a vector qvec within ixL
1080  !> Options to
1081  !> employ standard second order CD evaluations
1082  !> use Gauss theorem for non-Cartesian grids
1083  !> use Stokes theorem for non-Cartesian grids
1084  subroutine curlvector_trans(qvec,qvecc,ixI^L,ixO^L,curlvec,idim,idirmin,idirmin0,ndir0)
1086 
1087  integer, intent(in) :: ixI^L,ixO^L
1088  integer, intent(in) :: idim, ndir0, idirmin0
1089  integer, intent(inout) :: idirmin
1090  double precision, intent(in) :: qvec(ixi^s,1:ndir0),qvecc(ixi^s,1:ndir0)
1091  double precision, intent(inout) :: curlvec(ixi^s,idirmin0:3)
1092 
1093  integer :: ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L
1094  double precision :: invdx(1:ndim)
1095  double precision :: tmp(ixi^s),tmp2(ixi^s),xC(ixi^s),surface(ixi^s)
1096 
1097  ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k
1098  ! Curl can have components (idirmin:3)
1099  ! Determine exact value of idirmin while doing the loop.
1100 
1101  idirmin=4
1102  curlvec(ixo^s,idirmin0:3)=zero
1103 
1104  ! all non-Cartesian cases
1105  select case(coordinate)
1106  case(cartesian) ! Cartesian grids
1107  invdx=1.d0/dxlevel
1108  do idir=idirmin0,3
1109  do jdir=1,ndim; do kdir=1,ndir0
1110  if(lvc(idir,jdir,kdir)/=0)then
1111  if(jdir/=idim) then
1112  tmp(ixi^s)=qvec(ixi^s,kdir)
1113  hxo^l=ixo^l-kr(jdir,^d);
1114  jxo^l=ixo^l+kr(jdir,^d);
1115  else
1116  ! use two cell-center values for gradient at interface of the two cells
1117  ! because left and right neighbor interface values is unavailable at block boundary faces
1118  tmp(ixi^s)=qvecc(ixi^s,kdir)
1119  hxo^l=ixo^l;
1120  jxo^l=ixo^l+kr(jdir,^d);
1121  end if
1122  ! second order centered differencing
1123  tmp(ixo^s)=half*(tmp(jxo^s)-tmp(hxo^s))*invdx(jdir)
1124  !> \todo allow for 4th order CD evaluation here as well
1125  if(lvc(idir,jdir,kdir)==1)then
1126  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
1127  else
1128  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
1129  endif
1130  if(idir<idirmin)idirmin=idir
1131  endif
1132  enddo; enddo;
1133  end do
1134  case(cartesian_stretched) ! stretched Cartesian grids
1135  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
1136  if(lvc(idir,jdir,kdir)/=0)then
1137  select case(typecurl)
1138  case('central')
1139  tmp(ixi^s)=qvec(ixi^s,kdir)
1140  hxo^l=ixo^l-kr(jdir,^d);
1141  jxo^l=ixo^l+kr(jdir,^d);
1142  ! second order centered differencing
1143  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,jdir)-block%x(hxo^s,jdir))
1144  case('Gaussbased')
1145  hxo^l=ixo^l-kr(jdir,^d);
1146  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1147  jxc^l=ixc^l+kr(jdir,^d);
1148  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1149  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
1150  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
1151  case('Stokesbased')
1152  hxo^l=ixo^l-kr(jdir,^d);
1153  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1154  jxc^l=ixc^l+kr(jdir,^d);
1155  if(kdir<=ndim)then
1156  tmp(ixc^s)=block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1157  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
1158  else
1159  tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1160  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
1161  endif
1162  if(idir<=ndim)then
1163  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%surface(ixo^s,idir)
1164  else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
1165  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(block%ds(ixo^s,jdir)*block%ds(ixo^s,kdir))
1166  endif
1167  case default
1168  call mpistop('no such curl evaluator')
1169  end select
1170  if(lvc(idir,jdir,kdir)==1)then
1171  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1172  else
1173  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1174  endif
1175  if(idir<idirmin)idirmin=idir
1176  endif
1177  enddo; enddo; enddo;
1178  case(spherical) ! possibly stretched spherical grids
1179  select case(typecurl)
1180  case('central') ! ok for any dimensionality
1181  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
1182  if(lvc(idir,jdir,kdir)/=0)then
1183  tmp(ixi^s)=qvec(ixi^s,kdir)
1184  hxo^l=ixo^l-kr(jdir,^d);
1185  jxo^l=ixo^l+kr(jdir,^d);
1186  select case(jdir)
1187  case(1)
1188  tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1)
1189  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1))*block%x(ixo^s,1))
1190  {^nooned case(2)
1191  if(idir==1) tmp(ixi^s)=tmp(ixi^s)*dsin(block%x(ixi^s,2))
1192  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
1193  if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(block%x(ixo^s,2))
1194  }
1195  {^ifthreed case(3)
1196  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1)*dsin(block%x(ixo^s,2)))
1197  }
1198  end select
1199  if(lvc(idir,jdir,kdir)==1)then
1200  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1201  else
1202  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1203  endif
1204  if(idir<idirmin)idirmin=idir
1205  endif
1206  enddo; enddo; enddo;
1207  case('Gaussbased')
1208  do idir=idirmin0,3;
1209  do jdir=1,ndim; do kdir=1,ndir0
1210  if(lvc(idir,jdir,kdir)/=0)then
1211  hxo^l=ixo^l-kr(jdir,^d);
1212  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1213  jxc^l=ixc^l+kr(jdir,^d);
1214  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1215  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
1216  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
1217  if(lvc(idir,jdir,kdir)==1)then
1218  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1219  else
1220  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1221  endif
1222  if(idir<idirmin)idirmin=idir
1223  endif
1224  enddo; enddo;
1225  ! geometric terms
1226  if(idir==2.and.phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,phi_)/block%x(ixo^s,r_)
1227  {^nooned
1228  if(idir==phi_) curlvec(ixo^s,phi_)=curlvec(ixo^s,phi_)-qvec(ixo^s,2)/block%x(ixo^s,r_) &
1229  +qvec(ixo^s,r_)*dcos(block%x(ixo^s,2))/(block%x(ixo^s,r_)*dsin(block%x(ixo^s,2)))
1230  }
1231  enddo;
1232  case('Stokesbased')
1233  !if(ndim<3) call mpistop("Stokesbased for 3D spherical only")
1234  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
1235  if(lvc(idir,jdir,kdir)/=0)then
1236  select case(idir)
1237  case(1)
1238  if(jdir<kdir) then
1239  ! idir=1,jdir=2,kdir=3
1240  !! integral along 3rd dimension
1241  hxo^l=ixo^l-kr(jdir,^d);
1242  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1243  jxc^l=ixc^l+kr(jdir,^d);
1244  ! qvec(3) at cell interface along 2nd dimension
1245  if(stretched_dim(jdir) .and. stretch_uncentered) then
1246  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1247  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1248  else
1249  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1250  end if
1251  ! 2nd coordinate at cell interface along 2nd dimension
1252  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
1253  curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
1254  dsin(xc(hxo^s))*tmp(hxo^s))*block%dx(ixo^s,kdir)
1255  !! integral along 2nd dimension
1256  hxo^l=ixo^l-kr(kdir,^d);
1257  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1258  jxc^l=ixc^l+kr(kdir,^d);
1259  ! qvec(2) at cell interface along 3rd dimension
1260  if(stretched_dim(kdir) .and. stretch_uncentered) then
1261  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1262  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1263  else
1264  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1265  end if
1266  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir))&
1267  /block%surface(ixo^s,idir)*block%x(ixo^s,idir)
1268  end if
1269  case(2)
1270  if(jdir<kdir) then
1271  ! idir=2,jdir=1,kdir=3
1272  !! integral along 1st dimension
1273  hxo^l=ixo^l-kr(kdir,^d);
1274  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1275  jxc^l=ixc^l+kr(kdir,^d);
1276  ! qvec(1) at cell interface along 3rd dimension
1277  if(stretched_dim(kdir) .and. stretch_uncentered) then
1278  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1279  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1280  else
1281  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1282  end if
1283  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,1)
1284  !! integral along 3rd dimension
1285  hxo^l=ixo^l-kr(jdir,^d);
1286  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1287  jxc^l=ixc^l+kr(jdir,^d);
1288  ! qvec(3) at cell interface along 1st dimension
1289  if(stretched_dim(jdir) .and. stretch_uncentered) then
1290  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1291  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1292  else
1293  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1294  end if
1295  ! 1st coordinate at cell interface along 1st dimension
1296  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
1297  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
1298  dsin(block%x(ixo^s,idir))*block%dx(ixo^s,kdir))/block%surface(ixo^s,idir)
1299  end if
1300  case(3)
1301  if(jdir<kdir) then
1302  ! idir=3,jdir=1,kdir=2
1303  !! integral along 1st dimension
1304  hxo^l=ixo^l-kr(kdir,^d);
1305  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1306  jxc^l=ixc^l+kr(kdir,^d);
1307  ! qvec(1) at cell interface along 2nd dimension
1308  if(stretched_dim(kdir) .and. stretch_uncentered) then
1309  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1310  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1311  else
1312  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1313  end if
1314  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
1315  !! integral along 2nd dimension
1316  hxo^l=ixo^l-kr(jdir,^d);
1317  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1318  jxc^l=ixc^l+kr(jdir,^d);
1319  ! qvec(2) at cell interface along 1st dimension
1320  if(stretched_dim(jdir) .and. stretch_uncentered) then
1321  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1322  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1323  else
1324  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1325  end if
1326  ! 1st coordinate at cell interface along 1st dimension
1327  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
1328  if(ndim==3) then
1329  surface(ixo^s)=block%surface(ixo^s,idir)
1330  else
1331  surface(ixo^s)=block%x(ixo^s,jdir)*block%dx(ixo^s,kdir)*block%dx(ixo^s,jdir)
1332  end if
1333  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
1334  /surface(ixo^s)
1335  end if
1336  end select
1337  if(idir<idirmin)idirmin=idir
1338  endif
1339  enddo; enddo; enddo;
1340  case default
1341  call mpistop('no such curl evaluator')
1342  end select
1343  case(cylindrical) ! possibly stretched cylindrical grids
1344  select case(typecurl)
1345  case('central') ! works for any dimensionality, polar/cylindrical
1346  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
1347  if(lvc(idir,jdir,kdir)/=0)then
1348  tmp(ixi^s)=qvec(ixi^s,kdir)
1349  hxo^l=ixo^l-kr(jdir,^d);
1350  jxo^l=ixo^l+kr(jdir,^d);
1351  if(z_==3.or.z_==-1) then
1352  ! Case Polar_2D, Polar_2.5D or Polar_3D, i.e. R,phi,Z
1353  select case(jdir)
1354  case(1)
1355  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
1356  ! computes d(R V_phi)/dR or d V_Z/dR
1357  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
1358  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
1359  {^nooned case(2)
1360  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
1361  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
1362  }
1363  {^ifthreed case(3)
1364  ! handles d V_phi/dZ or d V_R/dZ
1365  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,3)-block%x(hxo^s,3))
1366  }
1367  end select
1368  end if
1369  if(phi_==3.or.phi_==-1) then
1370  ! Case Cylindrical_2D, Cylindrical_2.5D or Cylindrical_3D, i.e. R,Z,phi
1371  select case(jdir)
1372  case(1)
1373  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
1374  ! computes d(R V_phi)/dR or d V_Z/dR
1375  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
1376  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
1377  {^nooned case(2)
1378  ! handles d V_phi/dZ or d V_R/dZ
1379  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,2)-block%x(hxo^s,2))
1380  }
1381  {^ifthreed case(3)
1382  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
1383  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1))
1384  }
1385  end select
1386  end if
1387  if(lvc(idir,jdir,kdir)==1)then
1388  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1389  else
1390  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1391  endif
1392  if(idir<idirmin)idirmin=idir
1393  endif
1394  enddo; enddo; enddo;
1395  case('Gaussbased') ! works for any dimensionality, polar/cylindrical
1396  if(ndim<2) call mpistop("Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
1397  do idir=idirmin0,3;
1398  do jdir=1,ndim; do kdir=1,ndir0
1399  if(lvc(idir,jdir,kdir)/=0)then
1400  hxo^l=ixo^l-kr(jdir,^d);
1401  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1402  jxc^l=ixc^l+kr(jdir,^d);
1403  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1404  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
1405  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
1406  if(lvc(idir,jdir,kdir)==1)then
1407  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1408  else
1409  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1410  endif
1411  if(idir<idirmin)idirmin=idir
1412  endif
1413  enddo; enddo;
1414  ! geometric term from d e_R/d phi= e_phi for unit vectors e_R, e_phi
1415  ! but minus sign appears due to R,Z,phi ordering (?)
1416  ! note that in cylindrical 2D (RZ), phi_ is -1
1417  ! note that in polar 2D (Rphi), z_ is -1
1418  if((idir==phi_.or.(phi_==-1.and.idir==3)).and.z_>0) then
1419  ! cylindrical
1420  if( z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,z_)/block%x(ixo^s,r_)
1421  ! polar
1422  if(phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,z_)/block%x(ixo^s,r_)
1423  endif
1424  enddo;
1425  case('Stokesbased')
1426  if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only")
1427  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
1428  if(lvc(idir,jdir,kdir)/=0)then
1429  if(idir==r_) then
1430  if(jdir==phi_) then
1431  ! idir=r,jdir=phi,kdir=z
1432  !! integral along z dimension
1433  hxo^l=ixo^l-kr(jdir,^d);
1434  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1435  jxc^l=ixc^l+kr(jdir,^d);
1436  ! qvec(z) at cell interface along phi dimension
1437  if(stretched_dim(jdir) .and. stretch_uncentered) then
1438  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1439  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1440  else
1441  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1442  end if
1443  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,kdir)
1444  !! integral along phi dimension
1445  hxo^l=ixo^l-kr(kdir,^d);
1446  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1447  jxc^l=ixc^l+kr(kdir,^d);
1448  ! qvec(phi) at cell interface along z dimension
1449  if(stretched_dim(kdir) .and. stretch_uncentered) then
1450  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1451  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1452  else
1453  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1454  end if
1455  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%x(ixo^s,idir)*block%dx(ixo^s,jdir))&
1456  /block%surface(ixo^s,idir)
1457  end if
1458  else if(idir==phi_) then
1459  if(jdir<kdir) then
1460  ! idir=phi,jdir=r,kdir=z
1461  !! integral along r dimension
1462  hxo^l=ixo^l-kr(kdir,^d);
1463  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1464  jxc^l=ixc^l+kr(kdir,^d);
1465  ! qvec(r) at cell interface along z dimension
1466  if(stretched_dim(kdir) .and. stretch_uncentered) then
1467  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1468  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1469  else
1470  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1471  end if
1472  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,jdir)
1473  !! integral along z dimension
1474  hxo^l=ixo^l-kr(jdir,^d);
1475  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1476  jxc^l=ixc^l+kr(jdir,^d);
1477  ! qvec(z) at cell interface along r dimension
1478  if(stretched_dim(jdir) .and. stretch_uncentered) then
1479  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1480  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1481  else
1482  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1483  end if
1484  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,kdir))&
1485  /block%surface(ixo^s,idir)
1486  end if
1487  else ! idir==z_
1488  if(jdir<kdir) then
1489  ! idir=z,jdir=r,kdir=phi
1490  !! integral along r dimension
1491  hxo^l=ixo^l-kr(kdir,^d);
1492  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1493  jxc^l=ixc^l+kr(kdir,^d);
1494  ! qvec(r) at cell interface along phi dimension
1495  if(stretched_dim(kdir) .and. stretch_uncentered) then
1496  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
1497  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
1498  else
1499  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1500  end if
1501  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
1502  !! integral along phi dimension
1503  hxo^l=ixo^l-kr(jdir,^d);
1504  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
1505  jxc^l=ixc^l+kr(jdir,^d);
1506  ! qvec(phi) at cell interface along r dimension
1507  if(stretched_dim(jdir) .and. stretch_uncentered) then
1508  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
1509  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
1510  else
1511  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1512  end if
1513  ! r coordinate at cell interface along r dimension
1514  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
1515  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*block%dx(ixo^s,kdir))&
1516  /block%surface(ixo^s,idir)
1517  end if
1518  end if
1519  if(idir<idirmin)idirmin=idir
1520  endif
1521  enddo; enddo; enddo;
1522  case default
1523  call mpistop('no such curl evaluator')
1524  end select
1525  case default
1526  call mpistop('not possible to calculate curl')
1527  end select
1528 
1529  end subroutine curlvector_trans
1530 
1531 end module mod_geometry
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, parameter cartesian_stretched
Definition: mod_geometry.t:8
subroutine set_coordinate_system(geom)
Set the coordinate system to be used.
Definition: mod_geometry.t:17
procedure(set_surface), pointer usr_set_surface
subroutine set_pole
Definition: mod_geometry.t:95
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:310
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
Definition: mod_geometry.t:561
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
Definition: mod_geometry.t:354
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
integer coordinate
Definition: mod_geometry.t:6
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
logical, dimension(ndim) stretched_dim
True if a dimension is stretched.
integer, parameter limiter_ppm
Definition: mod_limiter.t:27
Module with slope/flux limiters.
Definition: mod_limiter.t:2
Module with all the methods that users can customize in AMRVAC.
character(len=std_len) geometry_name
logical, dimension(2, ndim) poleb
Indicates whether there is a pole at a boundary.
integer, parameter cartesian
Definition: mod_geometry.t:7
integer, parameter unitterm
Unit for standard output.
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:469
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision, dimension(:,:), allocatable dx
integer, dimension(:), allocatable, parameter d
subroutine curlvector_trans(qvec, qvecc, ixIL, ixOL, curlvec, idim, idirmin, idirmin0, ndir0)
Calculate idim transverse components of curl of a vector qvec within ixL Options to employ standard s...
integer mype
The rank of the current MPI task.
subroutine, public ppmlimitervar(ixIL, ixL, idims, q, qCT, qLC, qRC)
Definition: mod_ppm.t:12
Definition: mod_ppm.t:1
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
Definition: mod_limiter.t:128
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
subroutine get_surface_area(s, ixGL)
calculate area of surfaces of cells
Definition: mod_geometry.t:153
double precision, dimension(ndim) dxlevel
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stretch_uncentered
If true, adjust mod_geometry routines to account for grid stretching (but the flux computation will n...
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
Definition: mod_geometry.t:411
subroutine putgridgeo(igrid)
Deallocate geometry-related variables.
Definition: mod_geometry.t:143
logical, dimension(ndim) periodb
True for dimensions with periodic boundaries.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
Definition: mod_geometry.t:616
integer, parameter cartesian_expansion
Definition: mod_geometry.t:11
integer, parameter spherical
Definition: mod_geometry.t:10