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