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