MPI-AMRVAC  2.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 contains
7 
8  !> Set the coordinate system to be used
9  subroutine set_coordinate_system(geom)
11  character(len=*), intent(in) :: geom !< Name of the coordinate system
12 
13  select case (geom)
14  case ("Cartesian","Cartesian_1D","Cartesian_2D","Cartesian_3D")
15  ndir = ndim
16  typeaxial='slab'
17  case ("Cartesian_1.5D")
18  if (ndim /= 1) call mpistop("Geometry Cartesian_1.5D but ndim /= 1")
19  typeaxial='slab'
20  ndir = 2
21  case ("Cartesian_1.75D")
22  if (ndim /= 1) call mpistop("Geometry Cartesian_1.75D but ndim /= 1")
23  typeaxial='slab'
24  ndir = 3
25  case ("Cartesian_2.5D")
26  if (ndim /= 2) call mpistop("Geometry Cartesian_2.5D but ndim /= 2")
27  typeaxial='slab'
28  ndir = 3
29  case ("cylindrical","cylindrical_2D","cylindrical_3D")
30  ndir = ndim
31  r_ = 1
32  z_ = 2
33  if(ndir==3) phi_ = 3
34  typeaxial='cylindrical'
35  case ("cylindrical_2.5D")
36  if (ndim /= 2) call mpistop("Geometry cylindrical_2.5D but ndim /= 2")
37  ndir = 3
38  r_ = 1
39  z_ = 2
40  phi_ = 3
41  typeaxial='cylindrical'
42  case ("polar","polar_2D","polar_3D")
43  ndir = ndim
44  r_ = 1
45  phi_ = 2
46  if(ndir==3) z_ = 3
47  typeaxial='cylindrical'
48  case ("polar_1.5D")
49  if (ndim /= 1) call mpistop("Geometry polar_1.5D but ndim /= 1")
50  ndir = 2
51  r_ = 1
52  phi_ = 2
53  typeaxial='cylindrical'
54  case ("polar_2.5D")
55  if (ndim /= 2) call mpistop("Geometry polar_2.5D but ndim /= 2")
56  ndir = 3
57  r_ = 1
58  phi_ = 2
59  z_ = 3
60  typeaxial='cylindrical'
61  case ("spherical","spherical_2D","spherical_3D")
62  ndir = ndim
63  r_ = 1
64  if(ndir==3) phi_ = 3
65  z_ = -1
66  typeaxial='spherical'
67  case ("spherical_2.5D")
68  if (ndim /= 2) &
69  call mpistop("Geometry spherical_2.5D requires ndim == 2")
70  ndir = 3
71  r_ = 1
72  phi_ = 3
73  z_ = -1
74  typeaxial='spherical'
75  case default
76  call mpistop("Unknown geometry specified")
77  end select
78  end subroutine set_coordinate_system
79 
80  subroutine set_pole
82 
83  select case (typeaxial)
84  case ("spherical") {^ifthreed
85  ! For spherical grid, check whether phi-direction is periodic
86  if(periodb(ndim)) then
87  if(phi_/=3) call mpistop("phi_ should be 3 in 3D spherical coord!")
88  if(mod(ng3(1),2)/=0) &
89  call mpistop("Number of meshes in phi-direction should be even!")
90  if(abs(xprobmin2)<smalldouble) then
91  if(mype==0) write(unitterm,*) &
92  "Will apply pi-periodic conditions at northpole!"
93  poleb(1,2)=.true.
94  else
95  if(mype==0) write(unitterm,*) "There is no northpole!"
96  end if
97  if(abs(xprobmax2-dpi)<smalldouble) then
98  if(mype==0) write(unitterm,*) &
99  "Will apply pi-periodic conditions at southpole!"
100  poleb(2,2)=.true.
101  else
102  if(mype==0) write(unitterm,*) "There is no southpole!"
103  end if
104  end if}
105  case ("cylindrical")
106  {
107  if (^d == phi_ .and. periodb(^d)) then
108  if(mod(ng^d(1),2)/=0) then
109  call mpistop("Number of meshes in phi-direction should be even!")
110  end if
111 
112  if(abs(xprobmin1)<smalldouble) then
113  if (mype==0) then
114  write(unitterm,*) "Will apply pi-periodic conditions at r=0"
115  end if
116  poleb(1,1)=.true.
117  else
118  if (mype==0) then
119  write(unitterm,*) "There is no cylindrical axis!"
120  end if
121  end if
122  end if\}
123  end select
124 
125  end subroutine set_pole
126 
127  !> Deallocate geometry-related variables
128  subroutine putgridgeo(igrid)
130  integer, intent(in) :: igrid
131 
132  deallocate(ps(igrid)%surfaceC,ps(igrid)%surface,ps(igrid)%dvolume,ps(igrid)%dx,&
133  psc(igrid)%dx,ps(igrid)%ds,psc(igrid)%dvolume)
134 
135  end subroutine putgridgeo
136 
137  subroutine fillgeo(igrid,ixG^L)
139 
140  integer, intent(in) :: igrid, ixG^L
141 
142  integer :: ix^L, ixC^L
143  double precision :: x(ixg^s,ndim), drs(ixg^s), dx2(ixg^s), dx3(ixg^s)
144  !-----------------------------------------------------------------------------
145  ix^l=ixg^l^lsub1;
146 
147  select case (typeaxial)
148  case ("slabstretch")
149  drs(ixg^s)=ps(igrid)%dx(ixg^s,1)
150  {^nooned
151  dx2(ixg^s)=ps(igrid)%dx(ixg^s,2)}
152  {^ifthreed
153  dx3(ixg^s)=ps(igrid)%dx(ixg^s,3)}
154 
155  {^ifoned
156  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
157  ps(igrid)%surfaceC(ixc^s,1)=1.d0
158  ps(igrid)%surface(ixc^s,1) =1.d0
159  }
160  {^iftwod
161  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
162  ps(igrid)%surfaceC(ixc^s,1)=dx2(ixc^s)
163  ps(igrid)%surface(ixc^s,1) =dx2(ixc^s)
164  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
165  ps(igrid)%surfaceC(ixc^s,2)=drs(ixc^s)
166  ps(igrid)%surface(ixc^s,2)=drs(ixc^s)
167  }
168  {^ifthreed
169  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
170  ps(igrid)%surfaceC(ixc^s,1)= dx2(ixc^s)*dx3(ixc^s)
171  ps(igrid)%surface(ixc^s,1)=ps(igrid)%surfaceC(ixc^s,1)
172  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
173  ps(igrid)%surfaceC(ixc^s,2)= drs(ixc^s)*dx3(ixc^s)
174  ps(igrid)%surface(ixc^s,2)=ps(igrid)%surfaceC(ixc^s,2)
175  ixcmin^d=ixmin^d-kr(^d,3); ixcmax^d=ixmax^d;
176  ps(igrid)%surfaceC(ixc^s,3)= drs(ixc^s)*dx2(ixc^s)
177  ps(igrid)%surface(ixc^s,3)=ps(igrid)%surfaceC(ixc^s,3)
178  }
179 
180  case ("spherical")
181  x(ixg^s,1)=ps(igrid)%x(ixg^s,1)
182  {^nooned
183  x(ixg^s,2)=ps(igrid)%x(ixg^s,2)}
184 
185  drs(ixg^s)=ps(igrid)%dx(ixg^s,1)
186  {^nooned
187  dx2(ixg^s)=ps(igrid)%dx(ixg^s,2)}
188  {^ifthreed
189  dx3(ixg^s)=ps(igrid)%dx(ixg^s,3)}
190 
191  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
192 
193  ps(igrid)%surfaceC(ixc^s,1)=(x(ixc^s,1)+half*drs(ixc^s))**2 {^nooned &
194  *two*dsin(x(ixc^s,2))*dsin(half*dx2(ixc^s))}{^ifthreed*dx3(ixc^s)}
195 
196  {^nooned
197  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
198  ps(igrid)%surfaceC(ixc^s,2)=x(ixc^s,1)*drs(ixc^s)&
199  *dsin(x(ixc^s,2)+half*dx2(ixc^s))}{^ifthreed*dx3(ixc^s)}
200 
201  {^ifthreed
202  ixcmin^d=ixmin^d-kr(^d,3); ixcmax^d=ixmax^d;
203  ps(igrid)%surfaceC(ixc^s,3)=x(ixc^s,1)*drs(ixc^s)*dx2(ixc^s)}
204 
205  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
206  ps(igrid)%surface(ixc^s,1)=x(ixc^s,1)**2 {^nooned &
207  *two*dsin(x(ixc^s,2))*dsin(half*dx2(ixc^s))}{^ifthreed*dx3(ixc^s)}
208  {^nooned
209  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
210  ps(igrid)%surface(ixc^s,2)=x(ixc^s,1)*drs(ixc^s)&
211  *dsin(x(ixc^s,2))}{^ifthreed*dx3(ixc^s)}
212 
213  {^ifthreed
214  ixcmin^d=ixmin^d-kr(^d,3); ixcmax^d=ixmax^d;
215  ps(igrid)%surface(ixc^s,3)=x(ixc^s,1)*drs(ixc^s)*dx2(ixc^s)}
216 
217  case ("cylindrical")
218  x(ixg^s,1)=ps(igrid)%x(ixg^s,1)
219  drs(ixg^s)=ps(igrid)%dx(ixg^s,1)
220  {^nooned
221  dx2(ixg^s)=ps(igrid)%dx(ixg^s,2)}
222  {^ifthreed
223  dx3(ixg^s)=ps(igrid)%dx(ixg^s,3)}
224 
225  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
226  ps(igrid)%surfaceC(ixc^s,1)=dabs(x(ixc^s,1)+half*drs(ixc^s)){^de&*dx^de(ixc^s) }
227  {^nooned
228  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
229  if (z_==2) ps(igrid)%surfaceC(ixc^s,2)=x(ixc^s,1)*drs(ixc^s){^ifthreed*dx3(ixc^s)}
230  if (phi_==2) ps(igrid)%surfaceC(ixc^s,2)=drs(ixc^s){^ifthreed*dx3(ixc^s)}}
231  {^ifthreed
232  ixcmin^d=ixmin^d-kr(^d,3); ixcmax^d=ixmax^d;
233  if (z_==3) ps(igrid)%surfaceC(ixc^s,3)=x(ixc^s,1)*drs(ixc^s)*dx2(ixc^s)
234  if (phi_==3) ps(igrid)%surfaceC(ixc^s,3)=drs(ixc^s)*dx2(ixc^s)}
235 
236  ixcmin^d=ixmin^d-kr(^d,1); ixcmax^d=ixmax^d;
237  ps(igrid)%surface(ixc^s,1)=dabs(x(ixc^s,1)){^de&*dx^de(ixc^s) }
238  {^nooned
239  ixcmin^d=ixmin^d-kr(^d,2); ixcmax^d=ixmax^d;
240  if (z_==2) ps(igrid)%surface(ixc^s,2)=x(ixc^s,1)*drs(ixc^s){^ifthreed*dx3(ixc^s)}
241  if (phi_==2) ps(igrid)%surface(ixc^s,2)=drs(ixc^s){^ifthreed*dx3(ixc^s)}}
242  {^ifthreed
243  ixcmin^d=ixmin^d-kr(^d,3); ixcmax^d=ixmax^d;
244  if (z_==3) ps(igrid)%surface(ixc^s,3)=x(ixc^s,1)*drs(ixc^s)*dx2(ixc^s)
245  if (phi_==3) ps(igrid)%surface(ixc^s,3)=drs(ixc^s)*dx2(ixc^s)}
246 
247  case default
248  call mpistop("Sorry, typeaxial unknown")
249  end select
250 
251  end subroutine fillgeo
252 
253  !> Calculate gradient of a scalar q within ixL in direction idir
254  subroutine gradient(q,ixI^L,ixO^L,idir,gradq)
256 
257  integer, intent(in) :: ixI^L, ixO^L, idir
258  double precision, intent(in) :: q(ixi^s)
259  double precision, intent(inout) :: gradq(ixi^s)
260  double precision :: x(ixi^s,1:ndim)
261  integer :: jxO^L, hxO^L
262 
263  x(ixi^s,1:ndim)=block%x(ixi^s,1:ndim)
264 
265  hxo^l=ixo^l-kr(idir,^d);
266  jxo^l=ixo^l+kr(idir,^d);
267  if(slab) then
268  gradq(ixo^s)=half*(q(jxo^s)-q(hxo^s))/dxlevel(idir)
269  else
270  select case(typeaxial)
271  case('slabstretch')
272  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
273  case('spherical')
274  select case(idir)
275  case(1)
276  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,1)-x(hxo^s,1)))
277  {^nooned
278  case(2)
279  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
280  }
281  {^ifthreed
282  case(3)
283  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)))
284  }
285  end select
286  case('cylindrical')
287  if(idir==phi_) then
288  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,phi_)-x(hxo^s,phi_))*x(ixo^s,r_))
289  else
290  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
291  end if
292  case default
293  call mpistop('Unknown geometry')
294  end select
295  end if
296 
297  end subroutine gradient
298 
299  !> Calculate gradient of a scalar q within ixL in direction idir
300  !> first use limiter to go from cell center to edge
301  subroutine gradients(q,ixI^L,ixO^L,idir,gradq)
303  use mod_limiter
304  use mod_ppm
305 
306  integer, intent(in) :: ixI^L, ixO^L, idir
307  double precision, intent(in) :: q(ixi^s)
308  double precision, intent(inout) :: gradq(ixi^s)
309  double precision ,dimension(ixI^S) :: qC,qL,qR,dqC,ldq,rdq
310 
311  double precision :: x(ixi^s,1:ndim)
312  double precision :: invdx
313  integer :: hxO^L,ixC^L,jxC^L,gxC^L,hxC^L
314 
315  x(ixi^s,1:ndim)=block%x(ixi^s,1:ndim)
316 
317  invdx=1.d0/dxlevel(idir)
318  hxo^l=ixo^l-kr(idir,^d);
319  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
320  jxc^l=ixc^l+kr(idir,^d);
321  gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
322  hxc^l=gxc^l+kr(idir,^d);
323 
324  ! set the gradient limiter here
325  qr(gxc^s) = q(hxc^s)
326  ql(gxc^s) = q(gxc^s)
327  if (typegradlimiter/=limiter_ppm) then
328  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
329  call dwlimiter2(dqc,ixi^l,gxc^l,idir,typegradlimiter,ldw=ldq,rdw=rdq)
330  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
331  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
332  else
333  call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
334  endif
335 
336  if(slab) then
337  gradq(ixo^s)=half*(qr(ixo^s)-ql(hxo^s))*invdx
338  else
339  gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/block%dx(ixo^s,idir)
340  select case(typeaxial)
341  case('spherical')
342  select case(idir)
343  case(2)
344  gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
345  {^ifthreed
346  case(3)
347  gradq(ixo^s)=gradq(ixo^s)/(x(ixo^s,1)*dsin(x(ixo^s,2)))
348  }
349  end select
350  case('cylindrical')
351  if(idir==phi_) gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
352  end select
353  end if
354 
355  end subroutine gradients
356 
357  !> Calculate divergence of a vector qvec within ixL
358  subroutine divvector(qvec,ixI^L,ixO^L,divq, fourthorder)
360  integer, intent(in) :: ixI^L,ixO^L
361  double precision, intent(in) :: qvec(ixi^s,1:ndir)
362  double precision, intent(inout) :: divq(ixi^s)
363  logical, intent(in), optional :: fourthorder !< Default: false
364  logical :: use_4th_order
365  double precision :: qC(ixi^s), invdx(1:ndim)
366  integer :: jxO^L, hxO^L, ixC^L, jxC^L
367  integer :: idims, ix^L, gxO^L, kxO^L
368 
369  use_4th_order = .false.
370  if (present(fourthorder)) use_4th_order = fourthorder
371 
372  if (use_4th_order) then
373  if (.not. slab) &
374  call mpistop("divvector: 4th order only supported for slab geometry")
375  ! Fourth order, stencil width is two
376  ix^l=ixo^l^ladd2;
377  else
378  ! Second order, stencil width is one
379  ix^l=ixo^l^ladd1;
380  end if
381 
382  if (iximin^d>ixmin^d.or.iximax^d<ixmax^d|.or.) &
383  call mpistop("Error in divvector: Non-conforming input limits")
384 
385  invdx=1.d0/dxlevel
386  divq(ixo^s)=0.0d0
387 
388  if (slab) then
389  do idims=1,ndim
390  if (.not. use_4th_order) then
391  ! Use second order scheme
392  jxo^l=ixo^l+kr(idims,^d);
393  hxo^l=ixo^l-kr(idims,^d);
394  divq(ixo^s)=divq(ixo^s)+half*(qvec(jxo^s,idims) &
395  - qvec(hxo^s,idims))*invdx(idims)
396  else
397  ! Use fourth order scheme
398  kxo^l=ixo^l+2*kr(idims,^d);
399  jxo^l=ixo^l+kr(idims,^d);
400  hxo^l=ixo^l-kr(idims,^d);
401  gxo^l=ixo^l-2*kr(idims,^d);
402  divq(ixo^s)=divq(ixo^s)+&
403  (-qvec(kxo^s,idims) + 8.0d0 * qvec(jxo^s,idims) - 8.0d0 * &
404  qvec(hxo^s,idims) + qvec(gxo^s,idims))/(12.0d0 * dxlevel(idims))
405  end if
406  end do
407  else
408  do idims=1,ndim
409  hxo^l=ixo^l-kr(idims,^d);
410  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
411  jxc^l=ixc^l+kr(idims,^d);
412  if(stretched_dim(idims) .and. stretch_uncentered) then
413  ! linear interpolation at cell interface along stretched dimension
414  qc(ixc^s)=block%surfaceC(ixc^s,idims)*(qvec(ixc^s,idims)+0.5d0*block%dx(ixc^s,idims)*&
415  (qvec(jxc^s,idims)-qvec(ixc^s,idims))/(block%x(jxc^s,idims)-block%x(ixc^s,idims)))
416  else
417  qc(ixc^s)=block%surfaceC(ixc^s,idims)*half*(qvec(ixc^s,idims)+qvec(jxc^s,idims))
418  end if
419  divq(ixo^s)=divq(ixo^s)+qc(ixo^s)-qc(hxo^s)
420  end do
421  divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
422  end if
423 
424 
425  end subroutine divvector
426 
427  !> Calculate divergence of a vector qvec within ixL
428  !> using limited extrapolation to cell edges
429  subroutine divvectors(qvec,ixI^L,ixO^L,divq)
431  use mod_limiter
432  use mod_ppm
433 
434  integer, intent(in) :: ixI^L,ixO^L
435  double precision, intent(in) :: qvec(ixi^s,1:ndir)
436  double precision, intent(inout) :: divq(ixi^s)
437  double precision, dimension(ixI^S) :: qL,qR,dqC,ldq,rdq
438 
439  double precision :: invdx(1:ndim)
440  integer :: hxO^L,ixC^L,jxC^L,idims,ix^L,gxC^L,hxC^L
441 
442  ix^l=ixo^l^ladd2;
443 
444  if (iximin^d>ixmin^d.or.iximax^d<ixmax^d|.or.) &
445  call mpistop("Error in divvectorS: Non-conforming input limits")
446 
447  invdx=1.d0/dxlevel
448  divq(ixo^s)=zero
449  do idims=1,ndim
450  hxo^l=ixo^l-kr(idims,^d);
451  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
452  jxc^l=ixc^l+kr(idims,^d);
453  gxcmin^d=ixcmin^d-kr(idims,^d);gxcmax^d=jxcmax^d;
454  hxc^l=gxc^l+kr(idims,^d);
455  qr(gxc^s) = qvec(hxc^s,idims)
456  ql(gxc^s) = qvec(gxc^s,idims)
457  if(typegradlimiter/=limiter_ppm) then
458  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
459  call dwlimiter2(dqc,ixi^l,gxc^l,idims,typegradlimiter,ldw=ldq,rdw=rdq)
460  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
461  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
462  else
463  dqc(ixi^s)=qvec(ixi^s,idims)
464  call ppmlimitervar(ixi^l,ixo^l,idims,dqc,dqc,ql,qr)
465  endif
466 
467  if (slab) then
468  divq(ixo^s)=divq(ixo^s)+half*(qr(ixo^s)-ql(hxo^s))*invdx(idims)
469  else
470  qr(ixc^s)=block%surfaceC(ixc^s,idims)*qr(ixc^s)
471  ql(ixc^s)=block%surfaceC(ixc^s,idims)*ql(ixc^s)
472  divq(ixo^s)=divq(ixo^s)+qr(ixo^s)-ql(hxo^s)
473  end if
474  end do
475  if(.not.slab) divq(ixo^s)=divq(ixo^s)/block%dvolume(ixo^s)
476 
477  end subroutine divvectors
478 
479  !> Calculate curl of a vector qvec within ixL
480  !> Options to
481  !> employ standard second order CD evaluations
482  !> use Gauss theorem for non-Cartesian grids
483  !> use Stokes theorem for non-Cartesian grids
484  subroutine curlvector(qvec,ixI^L,ixO^L,curlvec,idirmin,idirmin0,ndir0)
486 
487  integer, intent(in) :: ixI^L,ixO^L
488  integer, intent(in) :: ndir0, idirmin0
489  integer, intent(inout) :: idirmin
490  double precision, intent(in) :: qvec(ixi^s,1:ndir0)
491  double precision, intent(inout) :: curlvec(ixi^s,idirmin0:3)
492 
493  integer :: ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L
494  double precision :: invdx(1:ndim)
495  double precision :: tmp(ixi^s),tmp2(ixi^s),xC(ixi^s),surface(ixi^s)
496  !-----------------------------------------------------------------------------
497 
498  ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k
499  ! Curl can have components (idirmin:3)
500  ! Determine exact value of idirmin while doing the loop.
501 
502  idirmin=4
503  curlvec(ixo^s,idirmin0:3)=zero
504 
505  if(slab) then ! Cartesian case
506  invdx=1.d0/dxlevel
507  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
508  if(lvc(idir,jdir,kdir)/=0)then
509  tmp(ixi^s)=qvec(ixi^s,kdir)
510  hxo^l=ixo^l-kr(jdir,^d);
511  jxo^l=ixo^l+kr(jdir,^d);
512  ! second order centered differencing
513  tmp(ixo^s)=half*(tmp(jxo^s)-tmp(hxo^s))*invdx(jdir)
514  !> \todo allow for 4th order CD evaluation here as well
515  if(lvc(idir,jdir,kdir)==1)then
516  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
517  else
518  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
519  endif
520  if(idir<idirmin)idirmin=idir
521  endif
522  enddo; enddo; enddo;
523  return
524  endif
525 
526  ! all non-Cartesian cases
527  select case(typeaxial)
528  case('slabstretch') ! stretched Cartesian grids
529  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
530  if(lvc(idir,jdir,kdir)/=0)then
531  select case(typecurl)
532  case('central')
533  tmp(ixi^s)=qvec(ixi^s,kdir)
534  hxo^l=ixo^l-kr(jdir,^d);
535  jxo^l=ixo^l+kr(jdir,^d);
536  ! second order centered differencing
537  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,jdir)-block%x(hxo^s,jdir))
538  case('Gaussbased')
539  hxo^l=ixo^l-kr(jdir,^d);
540  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
541  jxc^l=ixc^l+kr(jdir,^d);
542  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
543  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
544  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
545  case('Stokesbased')
546  hxo^l=ixo^l-kr(jdir,^d);
547  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
548  jxc^l=ixc^l+kr(jdir,^d);
549  if(kdir<=ndim)then
550  tmp(ixc^s)=block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
551  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
552  else
553  tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
554  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
555  endif
556  if(idir<=ndim)then
557  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%surface(ixo^s,idir)
558  else ! essentially for 2.5D case, idir=3 and jdir,kdir<=2
559  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(block%ds(ixo^s,jdir)*block%ds(ixo^s,kdir))
560  endif
561  case default
562  call mpistop('no such curl evaluator')
563  end select
564  if(lvc(idir,jdir,kdir)==1)then
565  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
566  else
567  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
568  endif
569  if(idir<idirmin)idirmin=idir
570  endif
571  enddo; enddo; enddo;
572  case('spherical') ! possibly stretched spherical grids
573  select case(typecurl)
574  case('central') ! ok for any dimensionality
575  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
576  if(lvc(idir,jdir,kdir)/=0)then
577  tmp(ixi^s)=qvec(ixi^s,kdir)
578  hxo^l=ixo^l-kr(jdir,^d);
579  jxo^l=ixo^l+kr(jdir,^d);
580  select case(jdir)
581  case(1)
582  tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1)
583  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,1)-block%x(hxo^s,1))*block%x(ixo^s,1))
584  {^nooned case(2)
585  if(idir==1) tmp(ixi^s)=tmp(ixi^s)*dsin(block%x(ixi^s,2))
586  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
587  if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(block%x(ixo^s,2))
588  }
589  {^ifthreed case(3)
590  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)))
591  }
592  end select
593  if(lvc(idir,jdir,kdir)==1)then
594  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
595  else
596  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
597  endif
598  if(idir<idirmin)idirmin=idir
599  endif
600  enddo; enddo; enddo;
601  case('Gaussbased')
602  do idir=idirmin0,3;
603  do jdir=1,ndim; do kdir=1,ndir0
604  if(lvc(idir,jdir,kdir)/=0)then
605  hxo^l=ixo^l-kr(jdir,^d);
606  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
607  jxc^l=ixc^l+kr(jdir,^d);
608  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
609  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
610  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
611  if(lvc(idir,jdir,kdir)==1)then
612  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
613  else
614  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
615  endif
616  if(idir<idirmin)idirmin=idir
617  endif
618  enddo; enddo;
619  ! geometric terms
620  if(idir==2.and.phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,phi_)/block%x(ixo^s,r_)
621  {^nooned
622  if(idir==phi_) curlvec(ixo^s,phi_)=curlvec(ixo^s,phi_)-qvec(ixo^s,2)/block%x(ixo^s,r_) &
623  +qvec(ixo^s,r_)*dcos(block%x(ixo^s,2))/(block%x(ixo^s,r_)*dsin(block%x(ixo^s,2)))
624  }
625  enddo;
626  case('Stokesbased')
627  !if(ndim<3) call mpistop("Stokesbased for 3D spherical only")
628  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
629  if(lvc(idir,jdir,kdir)/=0)then
630  select case(idir)
631  case(1)
632  if(jdir<kdir) then
633  ! idir=1,jdir=2,kdir=3
634  !! integral along 3rd dimension
635  hxo^l=ixo^l-kr(jdir,^d);
636  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
637  jxc^l=ixc^l+kr(jdir,^d);
638  ! qvec(3) at cell interface along 2nd dimension
639  if(stretched_dim(jdir) .and. stretch_uncentered) then
640  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
641  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
642  else
643  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
644  end if
645  ! 2nd coordinate at cell interface along 2nd dimension
646  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
647  curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
648  dsin(xc(hxo^s))*tmp(hxo^s))*block%dx(ixo^s,kdir)
649  !! integral along 2nd dimension
650  hxo^l=ixo^l-kr(kdir,^d);
651  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
652  jxc^l=ixc^l+kr(kdir,^d);
653  ! qvec(2) at cell interface along 3rd dimension
654  if(stretched_dim(kdir) .and. stretch_uncentered) then
655  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
656  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
657  else
658  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
659  end if
660  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir))&
661  /block%surface(ixo^s,idir)*block%x(ixo^s,idir)
662  end if
663  case(2)
664  if(jdir<kdir) then
665  ! idir=2,jdir=1,kdir=3
666  !! integral along 1st dimension
667  hxo^l=ixo^l-kr(kdir,^d);
668  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
669  jxc^l=ixc^l+kr(kdir,^d);
670  ! qvec(1) at cell interface along 3rd dimension
671  if(stretched_dim(kdir) .and. stretch_uncentered) then
672  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
673  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
674  else
675  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
676  end if
677  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,1)
678  !! integral along 3rd dimension
679  hxo^l=ixo^l-kr(jdir,^d);
680  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
681  jxc^l=ixc^l+kr(jdir,^d);
682  ! qvec(3) at cell interface along 1st dimension
683  if(stretched_dim(jdir) .and. stretch_uncentered) then
684  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
685  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
686  else
687  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
688  end if
689  ! 1st coordinate at cell interface along 1st dimension
690  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
691  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
692  dsin(block%x(ixo^s,idir))*block%dx(ixo^s,kdir))/block%surface(ixo^s,idir)
693  end if
694  case(3)
695  if(jdir<kdir) then
696  ! idir=3,jdir=1,kdir=2
697  !! integral along 1st dimension
698  hxo^l=ixo^l-kr(kdir,^d);
699  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
700  jxc^l=ixc^l+kr(kdir,^d);
701  ! qvec(1) at cell interface along 2nd dimension
702  if(stretched_dim(kdir) .and. stretch_uncentered) then
703  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
704  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
705  else
706  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
707  end if
708  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
709  !! integral along 2nd dimension
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  ! qvec(2) at cell interface along 1st dimension
714  if(stretched_dim(jdir) .and. stretch_uncentered) then
715  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
716  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
717  else
718  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
719  end if
720  ! 1st coordinate at cell interface along 1st dimension
721  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
722  if(ndim==3) then
723  surface(ixo^s)=block%surface(ixo^s,idir)
724  else
725  surface(ixo^s)=block%x(ixo^s,jdir)*block%dx(ixo^s,kdir)*block%dx(ixo^s,jdir)
726  end if
727  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))&
728  /surface(ixo^s)
729  end if
730  end select
731  if(idir<idirmin)idirmin=idir
732  endif
733  enddo; enddo; enddo;
734  case default
735  call mpistop('no such curl evaluator')
736  end select
737  case('cylindrical') ! possibly stretched cylindrical grids
738  select case(typecurl)
739  case('central') ! works for any dimensionality, polar/cylindrical
740  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
741  if(lvc(idir,jdir,kdir)/=0)then
742  tmp(ixi^s)=qvec(ixi^s,kdir)
743  hxo^l=ixo^l-kr(jdir,^d);
744  jxo^l=ixo^l+kr(jdir,^d);
745  if(z_==3.or.z_==-1) then
746  ! Case Polar_2D, Polar_2.5D or Polar_3D, i.e. R,phi,Z
747  select case(jdir)
748  case(1)
749  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
750  ! computes d(R V_phi)/dR or d V_Z/dR
751  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
752  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
753  {^nooned case(2)
754  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
755  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,2)-block%x(hxo^s,2))*block%x(ixo^s,1))
756  }
757  {^ifthreed case(3)
758  ! handles d V_phi/dZ or d V_R/dZ
759  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,3)-block%x(hxo^s,3))
760  }
761  end select
762  end if
763  if(phi_==3.or.phi_==-1) then
764  ! Case Cylindrical_2D, Cylindrical_2.5D or Cylindrical_3D, i.e. R,Z,phi
765  select case(jdir)
766  case(1)
767  if(idir==z_) tmp(ixi^s)=tmp(ixi^s)*block%x(ixi^s,1) ! R V_phi
768  ! computes d(R V_phi)/dR or d V_Z/dR
769  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,1)-block%x(hxo^s,1))
770  if(idir==z_) tmp2(ixo^s)=tmp2(ixo^s)/block%x(ixo^s,1) ! (1/R)*d(R V_phi)/dR
771  {^nooned case(2)
772  ! handles d V_phi/dZ or d V_R/dZ
773  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(block%x(jxo^s,2)-block%x(hxo^s,2))
774  }
775  {^ifthreed case(3)
776  ! handles (1/R)d V_Z/dphi or (1/R)d V_R/dphi
777  tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((block%x(jxo^s,3)-block%x(hxo^s,3))*block%x(ixo^s,1))
778  }
779  end select
780  end if
781  if(lvc(idir,jdir,kdir)==1)then
782  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
783  else
784  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
785  endif
786  if(idir<idirmin)idirmin=idir
787  endif
788  enddo; enddo; enddo;
789  case('Gaussbased') ! works for any dimensionality, polar/cylindrical
790  if(ndim<2) call mpistop("Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
791  do idir=idirmin0,3;
792  do jdir=1,ndim; do kdir=1,ndir0
793  if(lvc(idir,jdir,kdir)/=0)then
794  hxo^l=ixo^l-kr(jdir,^d);
795  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
796  jxc^l=ixc^l+kr(jdir,^d);
797  tmp(ixc^s)=block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
798  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir)))
799  tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/block%dvolume(ixo^s)
800  if(lvc(idir,jdir,kdir)==1)then
801  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
802  else
803  curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
804  endif
805  if(idir<idirmin)idirmin=idir
806  endif
807  enddo; enddo;
808  ! geometric term from d e_R/d phi= e_phi for unit vectors e_R, e_phi
809  ! but minus sign appears due to R,Z,phi ordering (?)
810  ! note that in cylindrical 2D (RZ), phi_ is -1
811  ! note that in polar 2D (Rphi), z_ is -1
812  if((idir==phi_.or.(phi_==-1.and.idir==3)).and.z_>0) then
813  ! cylindrical
814  if( z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,z_)/block%x(ixo^s,r_)
815  ! polar
816  if(phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,z_)/block%x(ixo^s,r_)
817  endif
818  enddo;
819  case('Stokesbased')
820  if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only")
821  do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0
822  if(lvc(idir,jdir,kdir)/=0)then
823  if(idir==r_) then
824  if(jdir==phi_) then
825  ! idir=r,jdir=phi,kdir=z
826  !! integral along z dimension
827  hxo^l=ixo^l-kr(jdir,^d);
828  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
829  jxc^l=ixc^l+kr(jdir,^d);
830  ! qvec(z) at cell interface along phi dimension
831  if(stretched_dim(jdir) .and. stretch_uncentered) then
832  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
833  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
834  else
835  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
836  end if
837  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,kdir)
838  !! integral along phi dimension
839  hxo^l=ixo^l-kr(kdir,^d);
840  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
841  jxc^l=ixc^l+kr(kdir,^d);
842  ! qvec(phi) at cell interface along z dimension
843  if(stretched_dim(kdir) .and. stretch_uncentered) then
844  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
845  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
846  else
847  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
848  end if
849  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%x(ixo^s,idir)*block%dx(ixo^s,jdir))&
850  /block%surface(ixo^s,idir)
851  end if
852  else if(idir==phi_) then
853  if(jdir<kdir) then
854  ! idir=phi,jdir=r,kdir=z
855  !! integral along r dimension
856  hxo^l=ixo^l-kr(kdir,^d);
857  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
858  jxc^l=ixc^l+kr(kdir,^d);
859  ! qvec(r) at cell interface along z dimension
860  if(stretched_dim(kdir) .and. stretch_uncentered) then
861  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
862  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
863  else
864  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
865  end if
866  curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*block%dx(ixo^s,jdir)
867  !! integral along z dimension
868  hxo^l=ixo^l-kr(jdir,^d);
869  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
870  jxc^l=ixc^l+kr(jdir,^d);
871  ! qvec(z) at cell interface along r dimension
872  if(stretched_dim(jdir) .and. stretch_uncentered) then
873  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
874  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
875  else
876  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
877  end if
878  curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,kdir))&
879  /block%surface(ixo^s,idir)
880  end if
881  else ! idir==z_
882  if(jdir<kdir) then
883  ! idir=z,jdir=r,kdir=phi
884  !! integral along r dimension
885  hxo^l=ixo^l-kr(kdir,^d);
886  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
887  jxc^l=ixc^l+kr(kdir,^d);
888  ! qvec(r) at cell interface along phi dimension
889  if(stretched_dim(kdir) .and. stretch_uncentered) then
890  tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*block%dx(ixc^s,kdir)*&
891  (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(block%x(jxc^s,kdir)-block%x(ixc^s,kdir))
892  else
893  tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
894  end if
895  curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*block%dx(ixo^s,jdir)
896  !! integral along phi dimension
897  hxo^l=ixo^l-kr(jdir,^d);
898  ixcmin^d=hxomin^d;ixcmax^d=ixomax^d;
899  jxc^l=ixc^l+kr(jdir,^d);
900  ! qvec(phi) at cell interface along r dimension
901  if(stretched_dim(jdir) .and. stretch_uncentered) then
902  tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*block%dx(ixc^s,jdir)*&
903  (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(block%x(jxc^s,jdir)-block%x(ixc^s,jdir))
904  else
905  tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
906  end if
907  ! r coordinate at cell interface along r dimension
908  xc(ixc^s)=block%x(ixc^s,jdir)+0.5d0*block%dx(ixc^s,jdir)
909  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))&
910  /block%surface(ixo^s,idir)
911  end if
912  end if
913  if(idir<idirmin)idirmin=idir
914  endif
915  enddo; enddo; enddo;
916  case default
917  call mpistop('no such curl evaluator')
918  end select
919  end select
920 
921  end subroutine curlvector
922 
923 end module mod_geometry
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine set_coordinate_system(geom)
Set the coordinate system to be used.
Definition: mod_geometry.t:10
subroutine set_pole
Definition: mod_geometry.t:81
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:255
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
Definition: mod_geometry.t:430
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine fillgeo(igrid, ixGL)
Definition: mod_geometry.t:138
character(len=std_len) typeaxial
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:22
Module with slope/flux limiters.
Definition: mod_limiter.t:2
logical, dimension(2, ndim) poleb
Indicates whether there is a pole at a boundary.
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:485
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 divvector(qvec, ixIL, ixOL, divq, fourthorder)
Calculate divergence of a vector qvec within ixL.
Definition: mod_geometry.t:359
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:160
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:302
logical slab
Cartesian geometry or not.
subroutine putgridgeo(igrid)
Deallocate geometry-related variables.
Definition: mod_geometry.t:129
logical, dimension(ndim) periodb
True for dimensions with periodic boundaries.
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:80