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