23 character(len=*),
intent(in) :: geom
29 case (
"Cartesian",
"Cartesian_1D",
"Cartesian_2D",
"Cartesian_3D")
32 case (
"Cartesian_1D_expansion")
33 if (
ndim /= 1)
call mpistop(
"Geometry Cartesian_1D_expansion but ndim /= 1")
36 case (
"Cartesian_1.5D")
37 if (
ndim /= 1)
call mpistop(
"Geometry Cartesian_1.5D but ndim /= 1")
40 case (
"Cartesian_1.75D")
41 if (
ndim /= 1)
call mpistop(
"Geometry Cartesian_1.75D but ndim /= 1")
44 case (
"Cartesian_2.5D")
45 if (
ndim /= 2)
call mpistop(
"Geometry Cartesian_2.5D but ndim /= 2")
48 case (
"cylindrical",
"cylindrical_2D",
"cylindrical_3D")
54 case (
"cylindrical_2.5D")
55 if (
ndim /= 2)
call mpistop(
"Geometry cylindrical_2.5D but ndim /= 2")
61 case (
"polar",
"polar_2D",
"polar_3D")
68 if (
ndim /= 1)
call mpistop(
"Geometry polar_1.5D but ndim /= 1")
74 if (
ndim /= 2)
call mpistop(
"Geometry polar_2.5D but ndim /= 2")
80 case (
"spherical",
"spherical_2D",
"spherical_3D")
86 case (
"spherical_2.5D")
88 call mpistop(
"Geometry spherical_2.5D requires ndim == 2")
95 call mpistop(
"Unknown geometry specified")
106 if(
phi_/=3)
call mpistop(
"phi_ should be 3 in 3D spherical coord!")
107 if(mod(ng3(1),2)/=0) &
108 call mpistop(
"Number of meshes in phi-direction should be even!")
109 if(abs(xprobmin2)<smalldouble)
then
111 "Will apply pi-periodic conditions at northpole!"
116 if(abs(xprobmax2-dpi)<smalldouble)
then
118 "Will apply pi-periodic conditions at southpole!"
126 if (^d == phi_ .and. periodb(^d))
then
127 if(mod(ng^d(1),2)/=0)
then
128 call mpistop(
"Number of meshes in phi-direction should be even!")
131 if(abs(xprobmin1)<smalldouble)
then
133 write(unitterm,*)
"Will apply pi-periodic conditions at r=0"
138 write(unitterm,*)
"There is no cylindrical axis!"
149 integer,
intent(in) :: igrid
151 deallocate(ps(igrid)%surfaceC,ps(igrid)%surface,ps(igrid)%dvolume,ps(igrid)%dx,&
152 psc(igrid)%dx,ps(igrid)%ds,psc(igrid)%dvolume)
162 integer,
intent(in) :: ixG^L
164 double precision :: x(ixG^S,ndim), xext(ixGmin1-1:ixGmax1,1), drs(ixG^S), drs_ext(ixGmin1-1:ixGmax1), dx2(ixG^S), dx3(ixG^S)
165 double precision :: exp_factor_ext(ixGmin1-1:ixGmax1),del_exp_factor_ext(ixGmin1-1:ixGmax1),exp_factor_primitive_ext(ixGmin1-1:ixGmax1)
166 double precision :: exp_factor(ixGmin1:ixGmax1),del_exp_factor(ixGmin1:ixGmax1),exp_factor_primitive(ixGmin1:ixGmax1)
171 drs(ixg^s)=s%dx(ixg^s,1)
172 x(ixg^s,1)=s%x(ixg^s,1)
175 call usr_set_surface(ixgmin1,ixgmax1,x,drs,exp_factor,del_exp_factor,exp_factor_primitive)
176 if (any(exp_factor <= zero))
call mpistop(
"The area must always be strictly positive!")
178 s%surface(ixg^s,1)=exp_factor(ixg^s)
179 xext(ixgmin1-1,1)=x(1,1)-half*drs(1)
180 xext(ixg^s,1)=x(ixg^s,1)+half*drs(ixg^s)
181 drs_ext(ixgmin1-1)=drs(1)
182 drs_ext(ixg^s)=drs(ixg^s)
184 s%surfaceC(ixgmin1-1:ixgmax1,1)=exp_factor_ext(ixgmin1-1:ixgmax1)
188 drs(ixg^s)=s%dx(ixg^s,1)
190 dx2(ixg^s)=s%dx(ixg^s,2)}
192 dx3(ixg^s)=s%dx(ixg^s,3)}
195 s%surfaceC(ixg^s,1)=1.d0
196 s%surface(ixg^s,1) =1.d0
199 s%surfaceC(ixg^s,1)=dx2(ixg^s)
200 s%surfaceC(ixg^s,2)=drs(ixg^s)
201 s%surface(ixg^s,1) =dx2(ixg^s)
202 s%surface(ixg^s,2)=drs(ixg^s)
205 s%surfaceC(ixg^s,1)= dx2(ixg^s)*dx3(ixg^s)
206 s%surfaceC(ixg^s,2)= drs(ixg^s)*dx3(ixg^s)
207 s%surfaceC(ixg^s,3)= drs(ixg^s)*dx2(ixg^s)
208 s%surface(ixg^s,1)=s%surfaceC(ixg^s,1)
209 s%surface(ixg^s,2)=s%surfaceC(ixg^s,2)
210 s%surface(ixg^s,3)=s%surfaceC(ixg^s,3)
212 {s%surfaceC(0^
d%ixG^s,^
d)=s%surfaceC(1^
d%ixG^s,^
d);\}
214 x(ixg^s,1)=s%x(ixg^s,1)
216 x(ixg^s,2)=s%x(ixg^s,2)
218 drs(ixg^s)=s%dx(ixg^s,1)
220 dx2(ixg^s)=s%dx(ixg^s,2)
223 dx3(ixg^s)=s%dx(ixg^s,3)
226 s%surfaceC(ixg^s,1)=(x(ixg^s,1)+half*drs(ixg^s))**2 {^nooned &
227 *two*dsin(x(ixg^s,2))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
230 s%surfaceC(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
231 *dsin(x(ixg^s,2)+half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
234 s%surfaceC(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
238 s%surfaceC(0,1)=dabs(x(1,1)-half*drs(1))**2
241 s%surfaceC(0,ixgmin2:ixgmax2,1)=(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
242 ixgmin2:ixgmax2))**2*two*dsin(x(1,ixgmin2:ixgmax2,2))*dsin(half*dx2(1,&
244 s%surfaceC(ixgmin1:ixgmax1,0,2)=x(ixgmin1:ixgmax1,1,&
245 1)*drs(ixgmin1:ixgmax1,1)*dsin(x(ixgmin1:ixgmax1,1,&
246 2)-half*dx2(ixgmin1:ixgmax1,1))
249 s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=(x(1,ixgmin2:ixgmax2,&
250 ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,&
251 ixgmin3:ixgmax3))**2*two*dsin(x(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3,&
252 2))*dsin(half*dx2(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx3(1,&
253 ixgmin2:ixgmax2,ixgmin3:ixgmax3)
254 s%surfaceC(ixgmin1:ixgmax1,0,ixgmin3:ixgmax3,2)=x(ixgmin1:ixgmax1,1,&
255 ixgmin3:ixgmax3,1)*drs(ixgmin1:ixgmax1,1,&
256 ixgmin3:ixgmax3)*dsin(x(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3,&
257 2)-half*dx2(ixgmin1:ixgmax1,1,ixgmin3:ixgmax3))*dx3(ixgmin1:ixgmax1,1,&
259 s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,0,3)=&
260 s%surfaceC(ixgmin1:ixgmax1,ixgmin2:ixgmax2,1,3)
263 s%surface(ixg^s,1)=x(ixg^s,1)**2 {^nooned &
264 *two*dsin(x(ixg^s,2))*dsin(half*dx2(ixg^s))}{^ifthreed*dx3(ixg^s)}
266 s%surface(ixg^s,2)=x(ixg^s,1)*drs(ixg^s)&
267 *dsin(x(ixg^s,2))}{^ifthreed*dx3(ixg^s)}
270 s%surface(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)}
273 x(ixg^s,1)=s%x(ixg^s,1)
274 drs(ixg^s)=s%dx(ixg^s,1)
276 dx2(ixg^s)=s%dx(ixg^s,2)}
278 dx3(ixg^s)=s%dx(ixg^s,3)}
280 s%surfaceC(ixg^s,1)=dabs(x(ixg^s,1)+half*drs(ixg^s)){^de&*
dx^de(ixg^s) }
282 if (
z_==2) s%surfaceC(ixg^s,2)=x(ixg^s,1)*drs(ixg^s){^ifthreed*dx3(ixg^s)}
283 if (
phi_==2) s%surfaceC(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}
286 if (
z_==3) s%surfaceC(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
287 if (
phi_==3) s%surfaceC(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)
290 s%surfaceC(0,1)=dabs(x(1,1)-half*drs(1))
293 s%surfaceC(0,ixgmin2:ixgmax2,1)=dabs(x(1,ixgmin2:ixgmax2,1)-half*drs(1,&
294 ixgmin2:ixgmax2))*dx2(1,ixgmin2:ixgmax2)
297 s%surfaceC(0,ixgmin2:ixgmax2,ixgmin3:ixgmax3,1)=dabs(x(1,ixgmin2:ixgmax2,&
298 ixgmin3:ixgmax3,1)-half*drs(1,ixgmin2:ixgmax2,ixgmin3:ixgmax3))*dx2(1,&
299 ixgmin2:ixgmax2,ixgmin3:ixgmax3)*dx3(1,ixgmin2:ixgmax2,&
302 {s%surfaceC(0^de%ixG^s,^de)=s%surfaceC(1^de%ixG^s,^de);\}
304 s%surface(ixg^s,1)=dabs(x(ixg^s,1)){^de&*
dx^de(ixg^s) }
306 if (
z_==2) s%surface(ixg^s,2)=x(ixg^s,1)*drs(ixg^s){^ifthreed*dx3(ixg^s)}
307 if (
phi_==2) s%surface(ixg^s,2)=drs(ixg^s){^ifthreed*dx3(ixg^s)}}
309 if (
z_==3) s%surface(ixg^s,3)=x(ixg^s,1)*drs(ixg^s)*dx2(ixg^s)
310 if (
phi_==3) s%surface(ixg^s,3)=drs(ixg^s)*dx2(ixg^s)}
313 call mpistop(
"Sorry, coordinate unknown")
322 integer,
intent(in) :: ixI^L, ixO^L, idir
323 double precision,
intent(in) :: q(ixI^S)
324 double precision,
intent(inout) :: gradq(ixI^S)
325 double precision :: x(ixI^S,1:ndim)
326 integer :: jxO^L, hxO^L
328 x(ixi^s,1:ndim)=
block%x(ixi^s,1:ndim)
330 hxo^l=ixo^l-
kr(idir,^
d);
331 jxo^l=ixo^l+
kr(idir,^
d);
334 gradq(ixo^s)=half*(q(jxo^s)-q(hxo^s))/
dxlevel(idir)
336 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
340 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,1)-x(hxo^s,1)))
343 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
347 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,3)-x(hxo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)))
352 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,
phi_)-x(hxo^s,
phi_))*x(ixo^s,
r_))
354 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
357 call mpistop(
'Unknown geometry')
363 subroutine gradientx(q,x,ixI^L,ixO^L,idir,gradq,fourth_order)
366 integer,
intent(in) :: ixI^L, ixO^L, idir
367 double precision,
intent(in) :: q(ixI^S), x(ixI^S,1:ndim)
368 double precision,
intent(inout) :: gradq(ixI^S)
369 logical,
intent(in) :: fourth_order
370 integer :: jxO^L, hxO^L, kxO^L
372 if(fourth_order)
then
375 if(iximin^
d>kxomin^
d.or.iximax^
d<kxomax^
d|.or.) &
376 call mpistop(
"Error in gradientx: Non-conforming input limits")
377 hxo^l=ixo^l-
kr(idir,^
d);
378 jxo^l=ixo^l+
kr(idir,^
d);
379 kxo^l=ixo^l+
kr(idir,^
d)*2;
383 jxo^l=ixo^l+
kr(idir,^
d);
386 if(fourth_order)
then
387 gradq(ixo^s)=(27.d0*(q(jxo^s)-q(ixo^s))-q(kxo^s)+q(hxo^s))/24.d0/
dxlevel(idir)
389 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/
dxlevel(idir)
392 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
396 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,1)-x(hxo^s,1)))
399 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,2)-x(hxo^s,2))*x(ixo^s,1))
403 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,3)-x(hxo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)))
408 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/((x(jxo^s,
phi_)-x(hxo^s,
phi_))*x(ixo^s,
r_))
410 gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
413 call mpistop(
'Unknown geometry')
425 integer,
intent(in) :: ixI^L, ixO^L, idir
426 double precision,
intent(in) :: q(ixI^S)
427 double precision,
intent(inout) :: gradq(ixI^S)
428 double precision ,
dimension(ixI^S) :: qC,qL,qR,dqC,ldq,rdq
430 double precision :: x(ixI^S,1:ndim)
431 double precision :: invdx
432 integer :: hxO^L,ixC^L,jxC^L,gxC^L,hxC^L
434 x(ixi^s,1:ndim)=
block%x(ixi^s,1:ndim)
437 hxo^l=ixo^l-
kr(idir,^
d);
438 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
439 jxc^l=ixc^l+
kr(idir,^
d);
440 gxcmin^
d=ixcmin^
d-
kr(idir,^
d);gxcmax^
d=jxcmax^
d;
441 hxc^l=gxc^l+
kr(idir,^
d);
447 dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
449 ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
450 qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
457 gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))*invdx
459 gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/
block%dx(ixo^s,idir)
461 gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/
block%dx(ixo^s,idir)
464 gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
467 gradq(ixo^s)=gradq(ixo^s)/(x(ixo^s,1)*dsin(x(ixo^s,2)))
471 gradq(ixo^s)=(qr(ixo^s)-ql(hxo^s))/
block%dx(ixo^s,idir)
472 if(idir==
phi_) gradq(ixo^s)=gradq(ixo^s)/x(ixo^s,1)
478 subroutine divvector(qvec,ixI^L,ixO^L,divq,fourthorder,sixthorder)
480 integer,
intent(in) :: ixI^L,ixO^L
481 double precision,
intent(in) :: qvec(ixI^S,1:ndir)
482 double precision,
intent(inout) :: divq(ixI^S)
483 logical,
intent(in),
optional :: fourthorder
484 logical,
intent(in),
optional :: sixthorder
485 logical :: use_4th_order
486 logical :: use_6th_order
487 double precision :: qC(ixI^S), invdx(1:ndim)
488 integer :: jxO^L, hxO^L, ixC^L, jxC^L
489 integer :: idims, ix^L, gxO^L, kxO^L
490 integer :: lxO^L, fxO^L
492 use_4th_order = .false.
493 use_6th_order = .false.
494 if (
present(fourthorder)) use_4th_order = fourthorder
495 if (
present(sixthorder)) use_6th_order = sixthorder
496 if(use_4th_order .and. use_6th_order) &
497 call mpistop(
"divvector: using 4th and 6th order at the same time")
499 if(use_4th_order)
then
501 call mpistop(
"divvector: 4th order only supported for slab geometry")
504 else if(use_6th_order)
then
507 call mpistop(
"divvector: 6th order only supported for slab geometry")
514 if (iximin^
d>ixmin^
d.or.iximax^
d<ixmax^
d|.or.) &
515 call mpistop(
"Error in divvector: Non-conforming input limits")
522 if(use_6th_order)
then
523 lxo^l=ixo^l+3*
kr(idims,^
d);
524 kxo^l=ixo^l+2*
kr(idims,^
d);
525 jxo^l=ixo^l+
kr(idims,^
d);
526 hxo^l=ixo^l-
kr(idims,^
d);
527 gxo^l=ixo^l-2*
kr(idims,^
d);
528 fxo^l=ixo^l-3*
kr(idims,^
d);
529 divq(ixo^s)=divq(ixo^s)+&
530 (-qvec(fxo^s,idims)+9.d0*qvec(gxo^s,idims)-45.d0*qvec(hxo^s,idims)&
531 +qvec(lxo^s,idims)-9.d0*qvec(kxo^s,idims)+45.d0*qvec(jxo^s,idims))&
533 else if(use_4th_order)
then
535 kxo^l=ixo^l+2*
kr(idims,^
d);
536 jxo^l=ixo^l+
kr(idims,^
d);
537 hxo^l=ixo^l-
kr(idims,^
d);
538 gxo^l=ixo^l-2*
kr(idims,^
d);
539 divq(ixo^s)=divq(ixo^s)+&
540 (-qvec(kxo^s,idims)+8.d0*qvec(jxo^s,idims)-8.d0*&
541 qvec(hxo^s,idims)+qvec(gxo^s,idims))/(12.d0*
dxlevel(idims))
544 jxo^l=ixo^l+
kr(idims,^
d);
545 hxo^l=ixo^l-
kr(idims,^
d);
546 divq(ixo^s)=divq(ixo^s)+half*(qvec(jxo^s,idims) &
547 - qvec(hxo^s,idims))*invdx(idims)
552 hxo^l=ixo^l-
kr(idims,^
d);
553 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
554 jxc^l=ixc^l+
kr(idims,^
d);
557 qc(ixc^s)=
block%surfaceC(ixc^s,idims)*(qvec(ixc^s,idims)+0.5d0*
block%dx(ixc^s,idims)*&
558 (qvec(jxc^s,idims)-qvec(ixc^s,idims))/(
block%x(jxc^s,idims)-
block%x(ixc^s,idims)))
560 qc(ixc^s)=
block%surfaceC(ixc^s,idims)*half*(qvec(ixc^s,idims)+qvec(jxc^s,idims))
562 divq(ixo^s)=divq(ixo^s)+qc(ixo^s)-qc(hxo^s)
564 divq(ixo^s)=divq(ixo^s)/
block%dvolume(ixo^s)
575 integer,
intent(in) :: ixI^L,ixO^L
576 double precision,
intent(in) :: qvec(ixI^S,1:ndir)
577 double precision,
intent(inout) :: divq(ixI^S)
578 double precision,
dimension(ixI^S) :: qL,qR,dqC,ldq,rdq
580 double precision :: invdx(1:ndim)
581 integer :: hxO^L,ixC^L,jxC^L,idims,ix^L,gxC^L,hxC^L
585 if (iximin^
d>ixmin^
d.or.iximax^
d<ixmax^
d|.or.) &
586 call mpistop(
"Error in divvectorS: Non-conforming input limits")
591 hxo^l=ixo^l-
kr(idims,^
d);
592 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
593 jxc^l=ixc^l+
kr(idims,^
d);
594 gxcmin^
d=ixcmin^
d-
kr(idims,^
d);gxcmax^
d=jxcmax^
d;
595 hxc^l=gxc^l+
kr(idims,^
d);
596 qr(gxc^s) = qvec(hxc^s,idims)
597 ql(gxc^s) = qvec(gxc^s,idims)
599 dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
601 ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
602 qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
604 dqc(ixi^s)=qvec(ixi^s,idims)
609 divq(ixo^s)=divq(ixo^s)+half*(qr(ixo^s)-ql(hxo^s))*invdx(idims)
611 qr(ixc^s)=
block%surfaceC(ixc^s,idims)*qr(ixc^s)
612 ql(ixc^s)=
block%surfaceC(ixc^s,idims)*ql(ixc^s)
613 divq(ixo^s)=divq(ixo^s)+qr(ixo^s)-ql(hxo^s)
625 subroutine curlvector(qvec,ixI^L,ixO^L,curlvec,idirmin,idirmin0,ndir0,fourthorder)
628 integer,
intent(in) :: ixI^L,ixO^L
629 integer,
intent(in) :: ndir0, idirmin0
630 integer,
intent(inout) :: idirmin
631 double precision,
intent(in) :: qvec(ixI^S,1:ndir0)
632 double precision,
intent(inout) :: curlvec(ixI^S,idirmin0:3)
633 logical,
intent(in),
optional :: fourthorder
635 integer :: ixA^L,ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L,kxO^L,gxO^L
636 double precision :: invdx(1:ndim)
637 double precision :: tmp(ixI^S),tmp2(ixI^S),xC(ixI^S),surface(ixI^S)
638 logical :: use_4th_order
644 use_4th_order = .false.
645 if (
present(fourthorder)) use_4th_order = fourthorder
647 if (use_4th_order)
then
649 call mpistop(
"curlvector: 4th order only supported for slab geometry")
657 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
658 call mpistop(
"Error in curlvector: Non-conforming input limits")
661 curlvec(ixo^s,idirmin0:3)=zero
667 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
668 if(
lvc(idir,jdir,kdir)/=0)
then
669 if (.not. use_4th_order)
then
671 jxo^l=ixo^l+
kr(jdir,^
d);
672 hxo^l=ixo^l-
kr(jdir,^
d);
673 tmp(ixo^s)=half*(qvec(jxo^s,kdir) &
674 - qvec(hxo^s,kdir))*invdx(jdir)
677 kxo^l=ixo^l+2*
kr(jdir,^
d);
678 jxo^l=ixo^l+
kr(jdir,^
d);
679 hxo^l=ixo^l-
kr(jdir,^
d);
680 gxo^l=ixo^l-2*
kr(jdir,^
d);
681 tmp(ixo^s)=(-qvec(kxo^s,kdir) + 8.0d0 * qvec(jxo^s,kdir) - 8.0d0 * &
682 qvec(hxo^s,kdir) + qvec(gxo^s,kdir))/(12.0d0 *
dxlevel(jdir))
684 if(
lvc(idir,jdir,kdir)==1)
then
685 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
687 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
689 if(idir<idirmin)idirmin=idir
693 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
694 if(
lvc(idir,jdir,kdir)/=0)
then
697 tmp(ixa^s)=qvec(ixa^s,kdir)
698 hxo^l=ixo^l-
kr(jdir,^
d);
699 jxo^l=ixo^l+
kr(jdir,^
d);
701 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,jdir)-
block%x(hxo^s,jdir))
703 hxo^l=ixo^l-
kr(jdir,^
d);
704 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
705 jxc^l=ixc^l+
kr(jdir,^
d);
706 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
707 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
708 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
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);
714 tmp(ixc^s)=
block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
715 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
717 tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
718 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
721 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%surface(ixo^s,idir)
723 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(
block%ds(ixo^s,jdir)*
block%ds(ixo^s,kdir))
726 call mpistop(
'no such curl evaluator')
728 if(
lvc(idir,jdir,kdir)==1)
then
729 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
731 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
733 if(idir<idirmin)idirmin=idir
739 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
740 if(
lvc(idir,jdir,kdir)/=0)
then
741 tmp(ixa^s)=qvec(ixa^s,kdir)
742 hxo^l=ixo^l-
kr(jdir,^
d);
743 jxo^l=ixo^l+
kr(jdir,^
d);
746 tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
747 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,1)-
block%x(hxo^s,1))*
block%x(ixo^s,1))
749 if(idir==1) tmp(ixa^s)=tmp(ixa^s)*dsin(
block%x(ixa^s,2))
750 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,2)-
block%x(hxo^s,2))*
block%x(ixo^s,1))
751 if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(
block%x(ixo^s,2))
754 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,3)-
block%x(hxo^s,3))*
block%x(ixo^s,1)*dsin(
block%x(ixo^s,2)))
757 if(
lvc(idir,jdir,kdir)==1)
then
758 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
760 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
762 if(idir<idirmin)idirmin=idir
767 do jdir=1,ndim;
do kdir=1,ndir0
768 if(
lvc(idir,jdir,kdir)/=0)
then
769 hxo^l=ixo^l-
kr(jdir,^
d);
770 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
771 jxc^l=ixc^l+
kr(jdir,^
d);
772 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
773 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
774 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
775 if(
lvc(idir,jdir,kdir)==1)
then
776 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
778 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
780 if(idir<idirmin)idirmin=idir
784 if(idir==2.and.
phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,
phi_)/
block%x(ixo^s,
r_)
792 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
793 if(
lvc(idir,jdir,kdir)/=0)
then
799 hxo^l=ixo^l-
kr(jdir,^
d);
800 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
801 jxc^l=ixc^l+
kr(jdir,^
d);
804 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
805 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
807 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
810 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
811 curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
812 dsin(xc(hxo^s))*tmp(hxo^s))*
block%dx(ixo^s,kdir)
814 hxo^l=ixo^l-
kr(kdir,^
d);
815 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
816 jxc^l=ixc^l+
kr(kdir,^
d);
819 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
820 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
822 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
824 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir))&
825 /
block%surface(ixo^s,idir)*
block%x(ixo^s,idir)
831 hxo^l=ixo^l-
kr(kdir,^
d);
832 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
833 jxc^l=ixc^l+
kr(kdir,^
d);
836 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
837 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
839 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
841 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,1)
843 hxo^l=ixo^l-
kr(jdir,^
d);
844 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
845 jxc^l=ixc^l+
kr(jdir,^
d);
848 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
849 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
851 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
854 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
855 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
856 dsin(
block%x(ixo^s,idir))*
block%dx(ixo^s,kdir))/
block%surface(ixo^s,idir)
862 hxo^l=ixo^l-
kr(kdir,^
d);
863 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
864 jxc^l=ixc^l+
kr(kdir,^
d);
867 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
868 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
870 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
872 curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir)
874 hxo^l=ixo^l-
kr(jdir,^
d);
875 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
876 jxc^l=ixc^l+
kr(jdir,^
d);
879 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
880 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
882 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
885 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
887 surface(ixo^s)=
block%surface(ixo^s,idir)
889 surface(ixo^s)=
block%x(ixo^s,jdir)*
block%dx(ixo^s,kdir)*
block%dx(ixo^s,jdir)
891 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*
block%dx(ixo^s,kdir))&
895 if(idir<idirmin)idirmin=idir
899 call mpistop(
'no such curl evaluator')
904 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
905 if(
lvc(idir,jdir,kdir)/=0)
then
906 tmp(ixa^s)=qvec(ixa^s,kdir)
907 hxo^l=ixo^l-
kr(jdir,^
d);
908 jxo^l=ixo^l+
kr(jdir,^
d);
909 if(
z_==3.or.
z_==-1)
then
913 if(idir==
z_) tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
915 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,1)-
block%x(hxo^s,1))
916 if(idir==
z_) tmp2(ixo^s)=tmp2(ixo^s)/
block%x(ixo^s,1)
919 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,2)-
block%x(hxo^s,2))*
block%x(ixo^s,1))
923 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,3)-
block%x(hxo^s,3))
931 if(idir==
z_) tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
933 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,1)-
block%x(hxo^s,1))
934 if(idir==
z_) tmp2(ixo^s)=tmp2(ixo^s)/
block%x(ixo^s,1)
937 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,2)-
block%x(hxo^s,2))
941 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,3)-
block%x(hxo^s,3))*
block%x(ixo^s,1))
945 if(
lvc(idir,jdir,kdir)==1)
then
946 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
948 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
950 if(idir<idirmin)idirmin=idir
954 if(ndim<2)
call mpistop(
"Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
956 do jdir=1,ndim;
do kdir=1,ndir0
957 if(
lvc(idir,jdir,kdir)/=0)
then
958 hxo^l=ixo^l-
kr(jdir,^
d);
959 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
960 jxc^l=ixc^l+
kr(jdir,^
d);
961 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
962 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
963 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
964 if(
lvc(idir,jdir,kdir)==1)
then
965 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
967 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
969 if(idir<idirmin)idirmin=idir
976 if((idir==
phi_.or.(
phi_==-1.and.idir==3)).and.
z_>0)
then
978 if(
z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,
z_)/
block%x(ixo^s,
r_)
980 if(
phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,
z_)/
block%x(ixo^s,
r_)
984 if(ndim<3)
call mpistop(
"Stokesbased for 3D cylindrical only")
985 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
986 if(
lvc(idir,jdir,kdir)/=0)
then
991 hxo^l=ixo^l-
kr(jdir,^
d);
992 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
993 jxc^l=ixc^l+
kr(jdir,^
d);
996 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
997 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
999 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1001 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,kdir)
1003 hxo^l=ixo^l-
kr(kdir,^
d);
1004 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1005 jxc^l=ixc^l+
kr(kdir,^
d);
1008 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1009 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1011 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1013 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%x(ixo^s,idir)*
block%dx(ixo^s,jdir))&
1014 /
block%surface(ixo^s,idir)
1016 else if(idir==
phi_)
then
1020 hxo^l=ixo^l-
kr(kdir,^
d);
1021 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1022 jxc^l=ixc^l+
kr(kdir,^
d);
1025 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1026 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1028 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1030 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,jdir)
1032 hxo^l=ixo^l-
kr(jdir,^
d);
1033 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1034 jxc^l=ixc^l+
kr(jdir,^
d);
1037 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1038 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1040 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1042 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,kdir))&
1043 /
block%surface(ixo^s,idir)
1049 hxo^l=ixo^l-
kr(kdir,^
d);
1050 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1051 jxc^l=ixc^l+
kr(kdir,^
d);
1054 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1055 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1057 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1059 curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir)
1061 hxo^l=ixo^l-
kr(jdir,^
d);
1062 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1063 jxc^l=ixc^l+
kr(jdir,^
d);
1066 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1067 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1069 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1072 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
1073 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*
block%dx(ixo^s,kdir))&
1074 /
block%surface(ixo^s,idir)
1077 if(idir<idirmin)idirmin=idir
1079 enddo; enddo; enddo;
1081 call mpistop(
'no such curl evaluator')
1084 call mpistop(
'not possible to calculate curl')
1097 integer,
intent(in) :: ixI^L,ixO^L
1098 integer,
intent(in) :: idim, ndir0, idirmin0
1099 integer,
intent(inout) :: idirmin
1100 double precision,
intent(in) :: qvec(ixI^S,1:ndir0),qvecc(ixI^S,1:ndir0)
1101 double precision,
intent(inout) :: curlvec(ixI^S,idirmin0:3)
1103 integer :: ixA^L,ixC^L,jxC^L,idir,jdir,kdir,hxO^L,jxO^L
1104 double precision :: invdx(1:ndim)
1105 double precision :: tmp(ixI^S),tmp2(ixI^S),xC(ixI^S),surface(ixI^S)
1112 curlvec(ixo^s,idirmin0:3)=zero
1121 do jdir=1,ndim;
do kdir=1,ndir0
1122 if(
lvc(idir,jdir,kdir)/=0)
then
1124 tmp(ixi^s)=qvec(ixi^s,kdir)
1125 hxo^l=ixo^l-
kr(jdir,^
d);
1126 jxo^l=ixo^l+
kr(jdir,^
d);
1130 tmp(ixi^s)=qvecc(ixi^s,kdir)
1132 jxo^l=ixo^l+
kr(jdir,^
d);
1135 tmp(ixo^s)=half*(tmp(jxo^s)-tmp(hxo^s))*invdx(jdir)
1136 if(
lvc(idir,jdir,kdir)==1)
then
1137 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp(ixo^s)
1139 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp(ixo^s)
1141 if(idir<idirmin)idirmin=idir
1146 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
1147 if(
lvc(idir,jdir,kdir)/=0)
then
1150 tmp(ixi^s)=qvec(ixi^s,kdir)
1151 hxo^l=ixo^l-
kr(jdir,^
d);
1152 jxo^l=ixo^l+
kr(jdir,^
d);
1154 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,jdir)-
block%x(hxo^s,jdir))
1156 hxo^l=ixo^l-
kr(jdir,^
d);
1157 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1158 jxc^l=ixc^l+
kr(jdir,^
d);
1159 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1160 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
1161 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
1163 hxo^l=ixo^l-
kr(jdir,^
d);
1164 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1165 jxc^l=ixc^l+
kr(jdir,^
d);
1167 tmp(ixc^s)=
block%ds(ixo^s,kdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1168 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
1170 tmp(ixc^s)=(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1171 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
1174 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%surface(ixo^s,idir)
1176 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/(
block%ds(ixo^s,jdir)*
block%ds(ixo^s,kdir))
1179 call mpistop(
'no such curl evaluator')
1181 if(
lvc(idir,jdir,kdir)==1)
then
1182 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1184 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1186 if(idir<idirmin)idirmin=idir
1188 enddo; enddo; enddo;
1192 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
1193 if(
lvc(idir,jdir,kdir)/=0)
then
1194 tmp(ixi^s)=qvec(ixi^s,kdir)
1195 hxo^l=ixo^l-
kr(jdir,^
d);
1196 jxo^l=ixo^l+
kr(jdir,^
d);
1199 tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
1200 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,1)-
block%x(hxo^s,1))*
block%x(ixo^s,1))
1202 if(idir==1) tmp(ixa^s)=tmp(ixa^s)*dsin(
block%x(ixa^s,2))
1203 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,2)-
block%x(hxo^s,2))*
block%x(ixo^s,1))
1204 if(idir==1) tmp2(ixo^s)=tmp2(ixo^s)/dsin(
block%x(ixo^s,2))
1207 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,3)-
block%x(hxo^s,3))*
block%x(ixo^s,1)*dsin(
block%x(ixo^s,2)))
1210 if(
lvc(idir,jdir,kdir)==1)
then
1211 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1213 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1215 if(idir<idirmin)idirmin=idir
1217 enddo; enddo; enddo;
1220 do jdir=1,ndim;
do kdir=1,ndir0
1221 if(
lvc(idir,jdir,kdir)/=0)
then
1222 hxo^l=ixo^l-
kr(jdir,^
d);
1223 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1224 jxc^l=ixc^l+
kr(jdir,^
d);
1225 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1226 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
1227 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
1228 if(
lvc(idir,jdir,kdir)==1)
then
1229 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1231 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1233 if(idir<idirmin)idirmin=idir
1237 if(idir==2.and.
phi_>0) curlvec(ixo^s,2)=curlvec(ixo^s,2)+qvec(ixo^s,
phi_)/
block%x(ixo^s,
r_)
1245 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
1246 if(
lvc(idir,jdir,kdir)/=0)
then
1252 hxo^l=ixo^l-
kr(jdir,^
d);
1253 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1254 jxc^l=ixc^l+
kr(jdir,^
d);
1257 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1258 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1260 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1263 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
1264 curlvec(ixo^s,idir)=(dsin(xc(ixo^s))*tmp(ixo^s)-&
1265 dsin(xc(hxo^s))*tmp(hxo^s))*
block%dx(ixo^s,kdir)
1267 hxo^l=ixo^l-
kr(kdir,^
d);
1268 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1269 jxc^l=ixc^l+
kr(kdir,^
d);
1272 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1273 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1275 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1277 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir))&
1278 /
block%surface(ixo^s,idir)*
block%x(ixo^s,idir)
1284 hxo^l=ixo^l-
kr(kdir,^
d);
1285 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1286 jxc^l=ixc^l+
kr(kdir,^
d);
1289 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1290 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1292 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1294 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,1)
1296 hxo^l=ixo^l-
kr(jdir,^
d);
1297 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1298 jxc^l=ixc^l+
kr(jdir,^
d);
1301 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1302 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1304 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1307 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
1308 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(hxo^s)*tmp(hxo^s)-xc(ixo^s)*tmp(ixo^s))*&
1309 dsin(
block%x(ixo^s,idir))*
block%dx(ixo^s,kdir))/
block%surface(ixo^s,idir)
1315 hxo^l=ixo^l-
kr(kdir,^
d);
1316 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1317 jxc^l=ixc^l+
kr(kdir,^
d);
1320 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1321 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1323 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1325 curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir)
1327 hxo^l=ixo^l-
kr(jdir,^
d);
1328 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1329 jxc^l=ixc^l+
kr(jdir,^
d);
1332 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1333 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1335 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1338 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
1340 surface(ixo^s)=
block%surface(ixo^s,idir)
1342 surface(ixo^s)=
block%x(ixo^s,jdir)*
block%dx(ixo^s,kdir)*
block%dx(ixo^s,jdir)
1344 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*
block%dx(ixo^s,kdir))&
1348 if(idir<idirmin)idirmin=idir
1350 enddo; enddo; enddo;
1352 call mpistop(
'no such curl evaluator')
1357 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
1358 if(
lvc(idir,jdir,kdir)/=0)
then
1359 tmp(ixi^s)=qvec(ixi^s,kdir)
1360 hxo^l=ixo^l-
kr(jdir,^
d);
1361 jxo^l=ixo^l+
kr(jdir,^
d);
1362 if(
z_==3.or.
z_==-1)
then
1366 if(idir==
z_) tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
1368 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,1)-
block%x(hxo^s,1))
1369 if(idir==
z_) tmp2(ixo^s)=tmp2(ixo^s)/
block%x(ixo^s,1)
1372 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,2)-
block%x(hxo^s,2))*
block%x(ixo^s,1))
1376 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,3)-
block%x(hxo^s,3))
1384 if(idir==
z_) tmp(ixa^s)=tmp(ixa^s)*
block%x(ixa^s,1)
1386 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,1)-
block%x(hxo^s,1))
1387 if(idir==
z_) tmp2(ixo^s)=tmp2(ixo^s)/
block%x(ixo^s,1)
1390 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/(
block%x(jxo^s,2)-
block%x(hxo^s,2))
1394 tmp2(ixo^s)=(tmp(jxo^s)-tmp(hxo^s))/((
block%x(jxo^s,3)-
block%x(hxo^s,3))*
block%x(ixo^s,1))
1398 if(
lvc(idir,jdir,kdir)==1)
then
1399 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1401 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1403 if(idir<idirmin)idirmin=idir
1405 enddo; enddo; enddo;
1407 if(ndim<2)
call mpistop(
"Gaussbased for 2D, 2.5D or 3D polar or cylindrical only")
1409 do jdir=1,ndim;
do kdir=1,ndir0
1410 if(
lvc(idir,jdir,kdir)/=0)
then
1411 hxo^l=ixo^l-
kr(jdir,^
d);
1412 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1413 jxc^l=ixc^l+
kr(jdir,^
d);
1414 tmp(ixc^s)=
block%surfaceC(ixc^s,jdir)*(qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1415 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir)))
1416 tmp2(ixo^s)=(tmp(ixo^s)-tmp(hxo^s))/
block%dvolume(ixo^s)
1417 if(
lvc(idir,jdir,kdir)==1)
then
1418 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+tmp2(ixo^s)
1420 curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-tmp2(ixo^s)
1422 if(idir<idirmin)idirmin=idir
1429 if((idir==
phi_.or.(
phi_==-1.and.idir==3)).and.
z_>0)
then
1431 if(
z_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)-qvec(ixo^s,
z_)/
block%x(ixo^s,
r_)
1433 if(
phi_==2) curlvec(ixo^s,idir)=curlvec(ixo^s,idir)+qvec(ixo^s,
z_)/
block%x(ixo^s,
r_)
1437 if(ndim<3)
call mpistop(
"Stokesbased for 3D cylindrical only")
1438 do idir=idirmin0,3;
do jdir=1,ndim;
do kdir=1,ndir0
1439 if(
lvc(idir,jdir,kdir)/=0)
then
1444 hxo^l=ixo^l-
kr(jdir,^
d);
1445 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1446 jxc^l=ixc^l+
kr(jdir,^
d);
1449 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1450 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1452 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1454 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,kdir)
1456 hxo^l=ixo^l-
kr(kdir,^
d);
1457 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1458 jxc^l=ixc^l+
kr(kdir,^
d);
1461 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1462 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1464 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1466 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%x(ixo^s,idir)*
block%dx(ixo^s,jdir))&
1467 /
block%surface(ixo^s,idir)
1469 else if(idir==
phi_)
then
1473 hxo^l=ixo^l-
kr(kdir,^
d);
1474 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1475 jxc^l=ixc^l+
kr(kdir,^
d);
1478 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1479 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1481 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1483 curlvec(ixo^s,idir)=(tmp(ixo^s)-tmp(hxo^s))*
block%dx(ixo^s,jdir)
1485 hxo^l=ixo^l-
kr(jdir,^
d);
1486 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1487 jxc^l=ixc^l+
kr(jdir,^
d);
1490 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1491 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1493 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1495 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,kdir))&
1496 /
block%surface(ixo^s,idir)
1502 hxo^l=ixo^l-
kr(kdir,^
d);
1503 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1504 jxc^l=ixc^l+
kr(kdir,^
d);
1507 tmp(ixc^s)=qvec(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,kdir)*&
1508 (qvec(jxc^s,jdir)-qvec(ixc^s,jdir))/(
block%x(jxc^s,kdir)-
block%x(ixc^s,kdir))
1510 tmp(ixc^s)=0.5d0*(qvec(ixc^s,jdir)+qvec(jxc^s,jdir))
1512 curlvec(ixo^s,idir)=(tmp(hxo^s)-tmp(ixo^s))*
block%dx(ixo^s,jdir)
1514 hxo^l=ixo^l-
kr(jdir,^
d);
1515 ixcmin^
d=hxomin^
d;ixcmax^
d=ixomax^
d;
1516 jxc^l=ixc^l+
kr(jdir,^
d);
1519 tmp(ixc^s)=qvec(ixc^s,kdir)+0.5d0*
block%dx(ixc^s,jdir)*&
1520 (qvec(jxc^s,kdir)-qvec(ixc^s,kdir))/(
block%x(jxc^s,jdir)-
block%x(ixc^s,jdir))
1522 tmp(ixc^s)=0.5d0*(qvec(ixc^s,kdir)+qvec(jxc^s,kdir))
1525 xc(ixc^s)=
block%x(ixc^s,jdir)+0.5d0*
block%dx(ixc^s,jdir)
1526 curlvec(ixo^s,idir)=(curlvec(ixo^s,idir)+(xc(ixo^s)*tmp(ixo^s)-xc(hxo^s)*tmp(hxo^s))*
block%dx(ixo^s,kdir))&
1527 /
block%surface(ixo^s,idir)
1530 if(idir<idirmin)idirmin=idir
1532 enddo; enddo; enddo;
1534 call mpistop(
'no such curl evaluator')
1537 call mpistop(
'not possible to calculate curl')
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Module with geometry-related routines (e.g., divergence, curl)
subroutine set_coordinate_system(geom)
Set the coordinate system to be used.
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
subroutine get_surface_area(s, ixGL)
calculate area of surfaces of cells
integer, parameter spherical
integer, parameter cartesian
integer, parameter stokesbased
integer, parameter cylindrical
integer, parameter gaussbased
integer, parameter cartesian_expansion
integer, parameter cartesian_stretched
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
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...
integer, parameter central
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...
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.
subroutine putgridgeo(igrid)
Deallocate geometry-related variables.
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
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.
integer, parameter limiter_ppm
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...
subroutine, public ppmlimitervar(ixIL, ixL, idims, q, qCT, qLC, qRC)
Module with all the methods that users can customize in AMRVAC.
procedure(set_surface), pointer usr_set_surface