!> Module with geometry-related routines (e.g., divergence, curl) module mod_geometry use mod_comm_lib, only: mpistop implicit none public integer :: coordinate=-1 integer, parameter :: Cartesian = 0 integer, parameter :: Cartesian_stretched= 1 integer, parameter :: cylindrical = 2 integer, parameter :: spherical = 3 integer, parameter :: Cartesian_expansion= 4 integer :: type_curl=0 integer, parameter :: central=1 integer, parameter :: Gaussbased=2 integer, parameter :: Stokesbased=3 contains !> Set the coordinate system from the GEOM compile-time define. !> !> Called once before usr_init, so a case does not have to call !> set_coordinate_system itself. Cases that do are still fine: the call is !> idempotent, and the consistency check in set_coordinate_system catches a !> name that disagrees with what the code was built for. !> !> Doing this centrally matters because set_coordinate_system is the only !> writer of coordinate, ndir, r_, phi_ and z_. Leaving coordinate at its !> -1 default is silently catastrophic rather than loud: read_par_files !> skips the conversion of angular domain bounds, and the select case in !> get_surface_area and in the cell-volume fill match nothing, so surfaceC !> and dvolume are simply never written. subroutine set_coordinate_system_from_config() call set_coordinate_system("Cartesian_3D") end subroutine set_coordinate_system_from_config !> Set the coordinate system to be used subroutine set_coordinate_system(geom) use mod_global_parameters character(len=*), intent(in) :: geom !< Name of the coordinate system ! Store the geometry name geometry_name = geom select case (geom) case ("Cartesian","Cartesian_1D","Cartesian_2D","Cartesian_3D") ndir = ndim coordinate=Cartesian case ("Cartesian_1D_expansion") if (ndim /= 1) call mpistop& ("Geometry Cartesian_1D_expansion but ndim /= 1") ndir = ndim coordinate=Cartesian_expansion case ("Cartesian_1.5D") if (ndim /= 1) call mpistop("Geometry Cartesian_1.5D but ndim /= 1") coordinate=Cartesian ndir = 2 case ("Cartesian_1.75D") if (ndim /= 1) call mpistop("Geometry Cartesian_1.75D but ndim /= 1") coordinate=Cartesian ndir = 3 case ("Cartesian_2.5D") if (ndim /= 2) call mpistop("Geometry Cartesian_2.5D but ndim /= 2") coordinate=Cartesian ndir = 3 case ("cylindrical","cylindrical_2D","cylindrical_3D") ndir = ndim r_ = 1 z_ = 2 if(ndir==3) phi_ = 3 coordinate=cylindrical case ("cylindrical_2.5D") if (ndim /= 2) call mpistop("Geometry cylindrical_2.5D but ndim /= 2") ndir = 3 r_ = 1 z_ = 2 phi_ = 3 coordinate=cylindrical case ("polar","polar_2D","polar_3D") ndir = ndim r_ = 1 phi_ = 2 if(ndir==3) z_ = 3 coordinate=cylindrical case ("polar_1.5D") if (ndim /= 1) call mpistop("Geometry polar_1.5D but ndim /= 1") ndir = 2 r_ = 1 phi_ = 2 coordinate=cylindrical case ("polar_2.5D") if (ndim /= 2) call mpistop("Geometry polar_2.5D but ndim /= 2") ndir = 3 r_ = 1 phi_ = 2 z_ = 3 coordinate=cylindrical case ("spherical","spherical_2D","spherical_3D") ndir = ndim r_ = 1 if(ndir==3) phi_ = 3 z_ = -1 coordinate=spherical case ("spherical_2.5D") if (ndim /= 2) call mpistop("Geometry spherical_2.5D requires ndim == 2") ndir = 3 r_ = 1 phi_ = 3 z_ = -1 coordinate=spherical ! The two logarithmically stretched systems are the same coordinate ! systems - identical ndir, r_, phi_, z_ and the same `coordinate` value, ! so every select case (coordinate) in the code keeps matching. What ! differs is only that the mesh is uniform in ln(r) rather than in r, which ! the LOG_RADIUS define carries. The names are kept distinct here because ! geometry_name goes into the snapshot header, and a log snapshot must not ! restart into a uniformly spaced build. case ("logSpherical","logSpherical_3D") ndir = ndim r_ = 1 if(ndir==3) phi_ = 3 z_ = -1 coordinate=spherical case ("logCylindrical","logCylindrical_3D") ndir = ndim r_ = 1 z_ = 2 if(ndir==3) phi_ = 3 coordinate=cylindrical case default call mpistop("Unknown geometry specified") end select ! The coordinate system is a compile-time choice in AGILE: the GEOM define ! (set from `geometry` in &meshlist) specializes the finite-volume kernels. ! Guard against a mod_usr that disagrees with what the code was built for. if (coordinate /= Cartesian) call mpistop("meshlist geometry is Cartesian, & &but mod_usr set a non-Cartesian coordinate system") ! coordinate alone cannot tell a log-stretched system from a uniform one - ! that is the point of the design - so the radial spacing is checked ! separately, on the name. if (geom(1:3) == "log") call mpistop("mod_usr set a logarithmically & &stretched coordinate system, but meshlist geometry is not & &logSpherical or logCylindrical") end subroutine set_coordinate_system !> The radial coordinate map, from the logical coordinate s to the physical !> radius r, and its inverse. !> !> These are the host-side definition of the map a LOG_RADIUS build uses, !> and the one place it is written down: !> !> log_r0 == 0 : r = exp(s) (s = ln(r)) !> log_r0 > 0 : r = sign(s)*r0*(exp(|s|) - 1) (s = ln(1 + r/r0)) !> !> The offset form is the one that reaches the axis. Its pivot is deliberate: !> s = 0 corresponds to r = 0 exactly, so the map is *odd*, and the mesh in !> the ghost layer beyond a cylindrical axis is the exact mirror of the mesh !> inside it. That is what the pole copy in getbc assumes when it fills ghost !> cell j from interior cell j of the block half a revolution away - the two !> have to occupy mirrored volumes, or the copied cell average belongs to a !> cell of a different size. The naive map r = exp(s) - r0 is *not* odd: !> its ghost cells shrink outward while the cells they mirror grow, and the !> first ghost cell is then misplaced by several percent of a cell width. !> !> The odd branch is only ever taken in that ghost layer; s >= 0 everywhere !> in the physical domain of an offset build. Device code that only ever !> sees the domain therefore uses the cheaper equivalent !> r = log_ra*exp(s) + log_rb instead of calling these. !> !> A case's mod_usr should use these rather than writing exp/log by hand: !> xprobmin1 and xprobmax1 hold the *logical* radial bounds in a LOG_RADIUS !> build, so r_of_s is how a case recovers a physical radius from them. pure elemental function r_of_s(s) result(r) use mod_global_parameters double precision, intent(in) :: s double precision :: r if (log_r0 > zero) then r = dsign(one,s)*log_r0*(dexp(dabs(s)) - one) else r = dexp(s) end if end function r_of_s !> Inverse of r_of_s. See there. pure elemental function s_of_r(r) result(s) use mod_global_parameters double precision, intent(in) :: r double precision :: s if (log_r0 > zero) then s = dsign(one,r)*dlog(one + dabs(r)/log_r0) else s = dlog(r) end if end function s_of_r subroutine set_pole use mod_global_parameters select case (coordinate) case (spherical) ! For spherical grid, check whether phi-direction is periodic if(periodB(ndim)) then if(phi_/=3) call mpistop("phi_ should be 3 in 3D spherical coord!") if(mod(ng3(1),2)/=0) call mpistop(& "Number of meshes in phi-direction should be even!") if(abs(xprobmin2) Validate a configured polar axis, once, on the host. !> !> Upstream MPI-AMRVAC defers most of this to the copy itself: pole_copy !> calls mpistop("Pole boundary condition should be symm or asymm") when it !> meets a typeboundary entry it cannot use. Here the equivalent branches !> live inside OpenACC kernels in getbc, where mpistop is unreachable and a !> bad entry would simply be treated as symmetric. So everything that could !> go wrong is checked here, before any of it is relied on. !> !> The checks also cover a mismatch that upstream cannot detect at all: where !> an axis is (poleB, derived from the domain bounds in set_pole) and where !> the user said one is (poleB_requested, recorded when read_par_files !> expanded typeboundary='pole') are decided independently, so a domain that !> reaches theta=0 without a 'pole' face - or a 'pole' face on a domain that !> does not reach the axis - is otherwise silently wrong. subroutine check_pole_setup use mod_global_parameters use mod_comm_lib, only: mpistop integer :: idim, iside, iB, iw double precision :: phi_len do idim=1,ndim do iside=1,2 if (poleB(iside,idim) .and. .not.poleB_requested(iside,idim)) then if (mype==0) write(unitterm,*) "dimension",idim," side",iside,& " reaches the axis but was not given typeboundary='pole'" call mpistop("a domain boundary lies on the axis but was not & &declared 'pole'") end if if (poleB_requested(iside,idim) .and. .not.poleB(iside,idim)) then if (mype==0) write(unitterm,*) "dimension",idim," side",iside,& " was given typeboundary='pole' but set_pole found no axis there" call mpistop("typeboundary='pole' on a face that is not on the & &axis; the phi direction has to be periodic, the number of level-1& & blocks across it even, and the face itself at theta=0, theta=pi & &or r=0") end if end do end do if (.not.any(poleB)) return if (stagger_grid) call mpistop("polar-axis treatment does not support & &stagger_grid") ! The pi-periodic neighbour is found by shifting the phi block index by ! ng(1)/2, which is half a revolution only if phi spans a full one. A ! half-turn domain with periodic boundaries passes every check upstream ! makes and then quietly shifts by the wrong angle. select case (phi_) case (1) phi_len = xprobmax1 - xprobmin1 case (2) phi_len = xprobmax2 - xprobmin2 case (3) phi_len = xprobmax3 - xprobmin3 case default call mpistop("a pole needs a phi direction") end select if (dabs(phi_len - two*dpi) > 1.0d-12*dpi) then if (mype==0) write(unitterm,*) "phi spans",phi_len," rather than 2*pi" call mpistop("a polar axis requires the phi direction to span a full & &turn (xprobmax-xprobmin = 1 in the 2*pi units of &meshlist)") end if ! getbc turns the per-variable entry on a pole face into a plus or a minus ! sign and has no third option, so anything else here would be read as ! symmetric. read_par_files writes only symm/asymm for a 'pole' face; this ! guards against a later edit that widens it. do idim=1,ndim do iside=1,2 if (.not.poleB(iside,idim)) cycle iB=2*idim-2+iside do iw=1,nwflux+nwaux if (typeboundary(iw,iB)/=bc_symm .and. typeboundary(iw,& iB)/=bc_asymm) then if (mype==0) write(unitterm,*) "variable",iw," on boundary",iB,& " has boundary type",typeboundary(iw,iB) call mpistop("every variable on a pole face must end up symm or & &asymm") end if end do end do end do end subroutine check_pole_setup !> Return cell-centre positions from the device to the host. !> !> fill_geometry_device builds a block's geometry on the GPU and leaves it !> there. Positions are the part host code asks for most often - the initial !> condition, the user hooks and the output all read ps(igrid)%x - but it !> asks rarely compared to how often blocks are (re)allocated, so the copy !> back does not belong in alloc_node: with AMR that would pay for it on !> every regrid, for blocks nothing on the host is going to look at. !> !> Call this instead immediately before host code that reads ps(igrid)%x, !> passing igrid when it is about to read a single block. Every such site is !> listed in CLAUDE.md; a missed one reads whatever the host copy held for !> the block that previously occupied that slot. subroutine sync_positions_host(igrid) use mod_global_parameters integer, intent(in), optional :: igrid integer :: iigrid if (present(igrid)) then else do iigrid = 1, igridstail end do end if end subroutine sync_positions_host !> Return a block's whole geometry - positions and metrics alike - from the !> device to the host, for every grid the host is about to read. !> !> This is the output-time counterpart of sync_positions_host: saveamrfile !> calls it in the same place it pulls back bg%w, which covers the log, the !> snapshot, the slices, the collapsed views and autoconvert, and the !> standalone convert path does the same before generate_plotfile. !> !> The fetch is unconditional and covers the whole grid in one pass, on !> purpose. Remembering "the host is already up to date" in a flag cannot be !> made honest, because the geometry is per block: any such flag would have !> to be cleared after refreshing whatever set of blocks happened to be !> listed at the time, and the next regrid or selectgrids can put a block !> that was not in that set straight back in front of host code. Doing every !> grid at once also keeps this to a few update directives per block rather !> than paying the per-call overhead once per reader. subroutine sync_geometry_host() use mod_global_parameters integer :: iigrid ! Index igrids(iigrid) directly rather than through a local igrid: a ! variable whose only uses are inside GPU directives looks dead to the ! Fortran front end, the assignment to it is dropped, and the update then ! silently targets the wrong block. sync_positions_host below is written ! the same way for the same reason. do iigrid = 1, igridstail end do end subroutine sync_geometry_host !> Calculate curl of a vector qvec within ixL !> Options to !> employ standard second order CD evaluations !> use Gauss theorem for non-Cartesian grids !> use Stokes theorem for non-Cartesian grids subroutine curlvector(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,& ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,curlvec,idirmin,idirmin0,& ndir0,fourthorder) use mod_global_parameters integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,& ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3 integer, intent(in) :: ndir0, idirmin0 integer, intent(inout) :: idirmin double precision, intent(in) :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3,1:ndir0) double precision, intent(inout) :: curlvec(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3,idirmin0:3) logical, intent(in), optional :: fourthorder !< Default: false integer :: ixAmin1,ixAmin2,ixAmin3,ixAmax1,ixAmax2,ixAmax3,& ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,jxCmin1,jxCmin2,jxCmin3,& jxCmax1,jxCmax2,jxCmax3,idir,jdir,kdir,hxOmin1,hxOmin2,hxOmin3,hxOmax1,& hxOmax2,hxOmax3,jxOmin1,jxOmin2,jxOmin3,jxOmax1,jxOmax2,jxOmax3,kxOmin1,& kxOmin2,kxOmin3,kxOmax1,kxOmax2,kxOmax3,gxOmin1,gxOmin2,gxOmin3,gxOmax1,& gxOmax2,gxOmax3 double precision :: invdx(1:ndim) double precision :: tmp(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),& tmp2(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),& xC(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),& surface(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3) logical :: use_4th_order ! Calculate curl within ixL: CurlV_i=eps_ijk*d_j V_k ! Curl can have components (idirmin:3) ! Determine exact value of idirmin while doing the loop. ! The Stokesbased discretisation divides by block%surface, the cell-centre ! face areas. Those are no longer built: no other routine reads them, and ! carrying a block-sized array per direction for 1:max_blocks purely for ! this branch is not worth the memory. Fail here rather than divide by an ! unassociated pointer. if (type_curl == Stokesbased) call mpistop(& "curlvector: typecurl='Stokesbased' needs block%surface, which is no & &longer built - see geo_t in mod_physicaldata") use_4th_order = .false. if (present(fourthorder)) use_4th_order = fourthorder if (use_4th_order) then if (.not. slab_uniform) call mpistop(& "curlvector: 4th order only supported for slab geometry") ! Fourth order, stencil width is two ixAmin1=ixOmin1-2;ixAmin2=ixOmin2-2;ixAmin3=ixOmin3-2;ixAmax1=ixOmax1+2 ixAmax2=ixOmax2+2;ixAmax3=ixOmax3+2; else ! Second order, stencil width is one ixAmin1=ixOmin1-1;ixAmin2=ixOmin2-1;ixAmin3=ixOmin3-1;ixAmax1=ixOmax1+1 ixAmax2=ixOmax2+1;ixAmax3=ixOmax3+1; end if if (ixImin1>ixAmin1.or.ixImax1ixAmin2.or.& ixImax2ixAmin3.or.ixImax30) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,2)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,2)+qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,phi_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_) if(idir==phi_) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,phi_)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,phi_)-qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,2)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_) +qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_)*dcos(block%x(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3,2))/(block%x(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3,& r_)*dsin(block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,2))) enddo; case(Stokesbased) !if(ndim<3) call mpistop("Stokesbased for 3D spherical only") do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0 if(lvc(idir,jdir,kdir)/=0)then select case(idir) case(1) if(jdir0) then ! cylindrical if( z_==2) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,idir)-qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,z_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_) ! polar if(phi_==2) curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,idir)=curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,idir)+qvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,z_)/block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_) endif enddo; case(Stokesbased) if(ndim<3) call mpistop("Stokesbased for 3D cylindrical only") do idir=idirmin0,3; do jdir=1,ndim; do kdir=1,ndir0 if(lvc(idir,jdir,kdir)/=0)then if(idir==r_) then if(jdir==phi_) then ! idir=r,jdir=phi,kdir=z !! integral along z dimension hxOmin1=ixOmin1-kr(jdir,1);hxOmin2=ixOmin2-kr(jdir,2) hxOmin3=ixOmin3-kr(jdir,3);hxOmax1=ixOmax1-kr(jdir,1) hxOmax2=ixOmax2-kr(jdir,2);hxOmax3=ixOmax3-kr(jdir,3); ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3 ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3; jxCmin1=ixCmin1+kr(jdir,1);jxCmin2=ixCmin2+kr(jdir,2) jxCmin3=ixCmin3+kr(jdir,3);jxCmax1=ixCmax1+kr(jdir,1) jxCmax2=ixCmax2+kr(jdir,2);jxCmax3=ixCmax3+kr(jdir,3); ! qvec(z) at cell interface along phi dimension if(stretched_dim(jdir) .and. stretch_uncentered) then tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3,kdir)+0.5d0*block%dx(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& jdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,kdir)-qvec(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& kdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,jdir)-block%x(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,jdir)) else tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& kdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,kdir)) end if curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,& idir)=(tmp(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)-tmp(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))*block%dx(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3,kdir) !! integral along phi dimension hxOmin1=ixOmin1-kr(kdir,1);hxOmin2=ixOmin2-kr(kdir,2) hxOmin3=ixOmin3-kr(kdir,3);hxOmax1=ixOmax1-kr(kdir,1) hxOmax2=ixOmax2-kr(kdir,2);hxOmax3=ixOmax3-kr(kdir,3); ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3 ixCmax1=ixOmax1;ixCmax2=ixOmax2;ixCmax3=ixOmax3; jxCmin1=ixCmin1+kr(kdir,1);jxCmin2=ixCmin2+kr(kdir,2) jxCmin3=ixCmin3+kr(kdir,3);jxCmax1=ixCmax1+kr(kdir,1) jxCmax2=ixCmax2+kr(kdir,2);jxCmax3=ixCmax3+kr(kdir,3); ! qvec(phi) at cell interface along z dimension if(stretched_dim(kdir) .and. stretch_uncentered) then tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3)=qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3,jdir)+0.5d0*block%dx(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& kdir)*(qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,jdir)-qvec(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& jdir))/(block%x(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,kdir)-block%x(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,kdir)) else tmp(ixCmin1:ixCmax1,ixCmin2:ixCmax2,& ixCmin3:ixCmax3)=0.5d0*(qvec(ixCmin1:ixCmax1,& ixCmin2:ixCmax2,ixCmin3:ixCmax3,& jdir)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,& jxCmin3:jxCmax3,jdir)) end if curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,& idir)=(curlvec(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,idir)+(tmp(hxOmin1:hxOmax1,& hxOmin2:hxOmax2,hxOmin3:hxOmax3)-tmp(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3))*block%x(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3,& idir)*block%dx(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,jdir))/block%surface(ixOmin1:ixOmax1,& ixOmin2:ixOmax2,ixOmin3:ixOmax3,idir) end if else if(idir==phi_) then if(jdir Calculate gradient of a scalar q within ixL in direction idir subroutine gradient(q,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,& ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,idir,gradq) use mod_global_parameters integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,& ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3, idir double precision, intent(in) :: q(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) double precision, intent(inout) :: gradq(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) integer :: jxOmin1,jxOmin2,jxOmin3,jxOmax1,jxOmax2,& jxOmax3, hxOmin1,hxOmin2,hxOmin3,hxOmax1,hxOmax2,hxOmax3 hxOmin1=ixOmin1-kr(idir,1);hxOmin2=ixOmin2-kr(idir,2) hxOmin3=ixOmin3-kr(idir,3);hxOmax1=ixOmax1-kr(idir,1) hxOmax2=ixOmax2-kr(idir,2);hxOmax3=ixOmax3-kr(idir,3); jxOmin1=ixOmin1+kr(idir,1);jxOmin2=ixOmin2+kr(idir,2) jxOmin3=ixOmin3+kr(idir,3);jxOmax1=ixOmax1+kr(idir,1) jxOmax2=ixOmax2+kr(idir,2);jxOmax3=ixOmax3+kr(idir,3); select case(coordinate) case(Cartesian) gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=half*(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/dxlevel(idir) case(Cartesian_stretched,Cartesian_expansion) gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,idir)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,idir)) case(spherical) select case(idir) case(1) gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,1)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,1))) case(2) gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,2)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,2))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,1)) case(3) gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,3)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,3))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,1)*dsin(block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,2))) end select case(cylindrical) if(idir==phi_) then gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/((block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,phi_)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,phi_))*block%x(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3,r_)) else gradq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,& ixOmin3:ixOmax3)=(q(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3)-q(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3))/(block%x(jxOmin1:jxOmax1,jxOmin2:jxOmax2,& jxOmin3:jxOmax3,idir)-block%x(hxOmin1:hxOmax1,hxOmin2:hxOmax2,& hxOmin3:hxOmax3,idir)) end if case default call mpistop('Unknown geometry') end select end subroutine gradient !> Calculate divergence of a vector qvec within ixL subroutine divvector(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,& ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,divq,fourthorder,& sixthorder) use mod_global_parameters integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,& ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3 double precision, intent(in) :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3,1:ndir) double precision, intent(inout) :: divq(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) logical, intent(in), optional :: fourthorder !< Default: false logical, intent(in), optional :: sixthorder !< Default: false logical :: use_4th_order logical :: use_6th_order double precision :: qC(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3), invdx(1:ndim) integer :: jxOmin1,jxOmin2,jxOmin3,jxOmax1,jxOmax2,& jxOmax3, hxOmin1,hxOmin2,hxOmin3,hxOmax1,hxOmax2,hxOmax3, ixCmin1,& ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3, jxCmin1,jxCmin2,jxCmin3,& jxCmax1,jxCmax2,jxCmax3 integer :: idims, ixmin1,ixmin2,ixmin3,ixmax1,& ixmax2,ixmax3, gxOmin1,gxOmin2,gxOmin3,gxOmax1,gxOmax2,gxOmax3, kxOmin1,& kxOmin2,kxOmin3,kxOmax1,kxOmax2,kxOmax3 integer :: lxOmin1,lxOmin2,lxOmin3,lxOmax1,lxOmax2,& lxOmax3, fxOmin1,fxOmin2,fxOmin3,fxOmax1,fxOmax2,fxOmax3 use_4th_order = .false. use_6th_order = .false. if (present(fourthorder)) use_4th_order = fourthorder if (present(sixthorder)) use_6th_order = sixthorder if(use_4th_order .and. use_6th_order) call & mpistop("divvector: using 4th and 6th order at the same time") if(use_4th_order) then if (.not. slab_uniform) call mpistop(& "divvector: 4th order only supported for slab geometry") ! Fourth order, stencil width is two ixmin1=ixOmin1-2;ixmin2=ixOmin2-2;ixmin3=ixOmin3-2;ixmax1=ixOmax1+2 ixmax2=ixOmax2+2;ixmax3=ixOmax3+2; else if(use_6th_order) then ! Sixth order, stencil width is three if (.not. slab_uniform) call mpistop(& "divvector: 6th order only supported for slab geometry") ixmin1=ixOmin1-3;ixmin2=ixOmin2-3;ixmin3=ixOmin3-3;ixmax1=ixOmax1+3 ixmax2=ixOmax2+3;ixmax3=ixOmax3+3; else ! Second order, stencil width is one ixmin1=ixOmin1-1;ixmin2=ixOmin2-1;ixmin3=ixOmin3-1;ixmax1=ixOmax1+1 ixmax2=ixOmax2+1;ixmax3=ixOmax3+1; end if if (ixImin1>ixmin1.or.ixImax1ixmin2.or.ixImax2ixmin3.or.ixImax3 Calculate divergence of a vector qvec within ixL !> using limited extrapolation to cell edges subroutine divvectorS(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,& ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,divq) use mod_global_parameters use mod_limiter integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,& ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3 double precision, intent(in) :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3,1:ndir) double precision, intent(inout) :: divq(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) double precision, dimension(ixImin1:ixImax1,ixImin2:ixImax2,& ixImin3:ixImax3) :: qL,qR,dqC,ldq,rdq double precision :: invdx(1:ndim) integer :: hxOmin1,hxOmin2,hxOmin3,hxOmax1,hxOmax2,hxOmax3,& ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,jxCmin1,jxCmin2,jxCmin3,& jxCmax1,jxCmax2,jxCmax3,idims,ixmin1,ixmin2,ixmin3,ixmax1,ixmax2,ixmax3,& gxCmin1,gxCmin2,gxCmin3,gxCmax1,gxCmax2,gxCmax3,hxCmin1,hxCmin2,hxCmin3,& hxCmax1,hxCmax2,hxCmax3 ixmin1=ixOmin1-2;ixmin2=ixOmin2-2;ixmin3=ixOmin3-2;ixmax1=ixOmax1+2 ixmax2=ixOmax2+2;ixmax3=ixOmax3+2; if (ixImin1>ixmin1.or.ixImax1ixmin2.or.ixImax2ixmin3.or.ixImax3