Available geometries

AGILE supports five coordinate settings: Cartesian (x, y, z), spherical (r, θ, φ), cylindrical (r, z, φ), and the two logarithmic-radius variants logSpherical and logCylindrical (the same two curvilinear systems with a radial mesh that is uniform in ln(r) rather than in r -- see Logarithmically stretched radius below). All are always 3D, so a curvilinear run is the full 3D system, never an axisymmetric or 1D/2D reduction.

The coordinate system is a compile-time choice: the finite-volume kernels are specialized for it, so that the (much cheaper) Cartesian kernel stays free of curvilinear bookkeeping. It is selected via geometry in &meshlist of the parameter file:

&meshlist
   geometry = 'spherical'   ! or 'cylindrical', 'logSpherical',
                            ! 'logCylindrical', or 'Cartesian' (the default)
/

logSpherical/logCylindrical are a spherical/cylindrical build (GEOM, and the runtime coordinate, stay 'spherical'/'cylindrical') plus a LOG_RADIUS flag; every geometry-specific branch stays correct untouched and only the radial spacing differs.

Conventions to know when setting up a curvilinear case:

• Angular domain bounds are given in units of 2π. So for spherical, xprobmin2 = 0.125d0, xprobmax2 = 0.375d0 places θ between π/4 and 3π/4. For cylindrical, only the φ bound (xprobmin3/xprobmax3) gets this treatment; the z bound (xprobmin2/xprobmax2) is a plain length, exactly like r.
• Axis order is (r, θ, φ) for spherical and (r, z, φ) for cylindrical -- in both cases φ is the third coordinate direction.
• Under logSpherical/logCylindrical, xprobmin1/xprobmax1 are physical radii in the parameter file and are converted to the logical ln-radius coordinate by read_par_files. A case's mod_usr.fpp that reads those bounds for a physical radius must map them back with r_of_s from mod_geometry (not by hand -- the map is only a plain exponential when log_r0 = 0).
• Cell positions ps(igrid)%x are the volume barycentre of the cell, not the midpoint of its faces -- so x ± dx/2 is not a cell face. See Cell positions are volume barycentres below.

Cell positions are volume barycentres

ps(igrid)%x -- which every case's mod_usr.fpp reads -- holds the volume barycentre of the cell. A cell average equals the point value at the barycentre to second order (the linear term of the Taylor expansion integrates to zero only about that point), so the barycentre is where the state, the initial condition and any position-dependent source term belong. Only the directions with a non-uniform volume weight are affected: r in both systems, and θ in spherical. φ in both, and z in cylindrical, carry uniform weight, so their barycentres are their midpoints. The offset is about 1% of a cell width radially on the test grids and grows to dθ/6 -- a sixth of a cell, pushed away from the axis -- for the first polar cell off θ = 0.

fill_geometry_device is written face-first: the radial-cell macro produces the two faces, the midpoint and the extent, and every area and volume is then written in terms of those so the expressions stay algebraically exact for any face pair. Three places in the tree that assumed x ± dx/2 was a face -- the face coordinates handed to the Riemann solver in mod_finite_volume, the VTU corner grid in calc_x, and fill_geometry_device itself -- rebuild the face from rnode instead. bgeo%ds deliberately stays on the cell midpoint: it is a geometric extent that setdt turns into the CFL length, and following the barycentre there would quietly relax the timestep.

Why curvilinear coordinates need extra source terms

AGILE discretizes the generic conservation law

$$\begin{equation} \partial_t \mathbf{U} + \nabla \cdot \mathbf{F}(\mathbf{U}) = \mathbf{S}(\mathbf{U},\nabla \mathbf{U},\mathbf{x},t) \end{equation}$$

(see Physics modules and equations) with a finite-volume update that divides the net flux through a cell's faces by the cell volume. In spherical or cylindrical coordinates the face areas and cell volume depend on position (e.g. a spherical shell's outer face is larger than its inner face), so this ratio is what actually implements the curvilinear divergence for scalar quantities -- no extra term is needed for the mass, energy, or (for ffhd) the field-aligned quantities, since these are advected as genuine scalars and Gauss's theorem then guarantees the finite volume update alone reproduces $\nabla \cdot \mathbf{F}$ correctly in any coordinate system.

Momentum (and, for mhd, the magnetic field) is different: it is a vector quantity, and its flux is a tensor built from a position-dependent (r, θ, φ) or (r, z, φ) basis. Writing the vector Laplacian/divergence out in curvilinear coordinates leaves over extra terms beyond the flux-divergence that the finite-volume update captures -- physically, these account for how the basis vectors themselves change from one face of a cell to the other. AGILE adds these leftover terms explicitly as a per-cell geometric source, in the addsource_geometry routine of each physics module (src/{hd,mhd,srhd,ffhd}/mod_*_templates.fpp); a physics module that does not implement one aborts at startup rather than silently omitting these terms. The explicit expressions below are exactly what these routines evaluate (rewritten into physical notation matching Physics modules and equations; see the source files for the exact discretization, in particular how the isotropic-pressure-like prefactor is evaluated from the actual discrete face areas rather than plugging in the continuous coordinate factor directly, so that a uniform pressure produces no net force to the same precision the code otherwise resolves -- for cylindrical this recovers the continuous prefactor exactly, since the radial face area of a cylindrical cell is linear in r). Because those discrete factors are the exact volume averages of the continuous 2/r, cot(θ)/r (spherical) or 1/r (cylindrical) prefactors for any face pair, no curvature term depends on where in the cell the stored position sits, and addsource_geometry takes no position argument -- which is what keeps it correct now that ps%x is the volume barycentre rather than the face midpoint.

Geometric source terms

The equations below give the additional source term $\mathbf{S}_{\mathrm{geom}}$ that a curvilinear run adds to the corresponding equation from Physics modules and equations; everywhere else the equations are unchanged. ffhd needs none, for any geometry (see above), so it is omitted below beyond noting $\mathbf{S}_{\mathrm{geom}} = 0$.

Spherical (r, θ, φ)

Hydrodynamics: hd

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{2p+\rho\left(v_\theta^2+v_\varphi^2\right)}{r} \,,\\ S_{\mathrm{geom}}[m_\theta] &=& \frac{p\cot\theta+\rho\left(v_\varphi^2\cot\theta-v_r v_\theta\right)}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& -\frac{\rho v_\varphi\left(v_r+v_\theta\cot\theta\right)}{r} \,. \end{eqnarray}$$

Magnetohydrodynamics: mhd

With the total (thermal + magnetic) pressure $p_{\mathrm{tot}} = p + \frac{1}{2}B^2$:

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{2p_{\mathrm{tot}}+\rho\left(v_\theta^2+v_\varphi^2\right)-\left(B_\theta^2+B_\varphi^2\right)}{r} \,,\\ S_{\mathrm{geom}}[m_\theta] &=& \frac{p_{\mathrm{tot}}\cot\theta+\cot\theta\left(\rho v_\varphi^2-B_\varphi^2\right)-\rho v_r v_\theta+B_r B_\theta}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& \frac{-\rho v_\varphi\left(v_r+v_\theta\cot\theta\right)+B_\varphi\left(B_r+B_\theta\cot\theta\right)}{r} \,,\\ S_{\mathrm{geom}}[B_r] &=& \frac{2\psi}{r} \,,\\ S_{\mathrm{geom}}[B_\theta] &=& \frac{v_r B_\theta-v_\theta B_r+\psi\cot\theta}{r} \,,\\ S_{\mathrm{geom}}[B_\varphi] &=& \frac{v_r B_\varphi-v_\varphi B_r+\cot\theta\left(v_\theta B_\varphi-v_\varphi B_\theta\right)}{r} \,, \end{eqnarray}$$

where ψ is the GLM field (see Divergence B source treatments); ψ itself has no geometric source.

Special relativistic hydrodynamics: srhd

Structurally identical to the hd case, with the momentum density $\mathbf{m} = \xi\mathbf{v}$ (rather than $\rho\mathbf{v}$) and using the spatial four-velocity components $\mathbf{u}\equiv\Gamma\mathbf{v}$ (the primitive momentum-slot variables; $\mathbf{v}=\mathbf{u}/\Gamma$), so that $\xi v_iv_j = (\xi/\Gamma^2)\,u_iu_j$:

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{2p+(\xi/\Gamma^2)\left(u_\theta^2+u_\varphi^2\right)}{r} \,,\\ S_{\mathrm{geom}}[m_\theta] &=& \frac{p\cot\theta+(\xi/\Gamma^2)\left(u_\varphi^2\cot\theta-u_r u_\theta\right)}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& -\frac{(\xi/\Gamma^2)\,u_\varphi\left(u_r+u_\theta\cot\theta\right)}{r} \,, \end{eqnarray}$$

with $\xi=\Gamma^2\rho h$ as in Physics modules and equations.

Cylindrical (r, z, φ)

Only the $m_r$/$m_\varphi$ (and, for mhd, $B_r$/$B_\varphi$) components pick up a geometric source; $m_z$/$B_z$ do not, since a cylindrical volume element's z-extent does not depend on r.

Hydrodynamics: hd

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{p+\rho v_\varphi^2}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& -\frac{\rho v_\varphi v_r}{r} \,. \end{eqnarray}$$

Magnetohydrodynamics: mhd

With $p_{\mathrm{tot}} = p + \frac{1}{2}B^2$:

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{p_{\mathrm{tot}}-B_\varphi^2+\rho v_\varphi^2}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& \frac{-\rho v_\varphi v_r+B_\varphi B_r}{r} \,,\\ S_{\mathrm{geom}}[B_r] &=& \frac{\psi}{r} \,,\\ S_{\mathrm{geom}}[B_\varphi] &=& \frac{B_\varphi v_r-B_r v_\varphi}{r} \,. \end{eqnarray}$$

Special relativistic hydrodynamics: srhd

As for spherical, with $\mathbf{u}\equiv\Gamma\mathbf{v}$ and $\xi=\Gamma^2\rho h$:

$$\begin{eqnarray} S_{\mathrm{geom}}[m_r] &=& \frac{p+(\xi/\Gamma^2)\,u_\varphi^2}{r} \,,\\ S_{\mathrm{geom}}[m_\varphi] &=& -\frac{(\xi/\Gamma^2)\,u_\varphi u_r}{r} \,. \end{eqnarray}$$

Logarithmically stretched radius

geometry = 'logSpherical' and geometry = 'logCylindrical' place the radial mesh uniformly in ln(r + r0) rather than in r. The offset r0 is the log_r0 entry of &meshlist and defaults to zero, the plain ln(r): the cell width then grows in proportion to r, so the relative resolution dr/r is constant -- the one stretching astrophysics usually wants, and the only one supported.

The design principle is that taking the logical radial coordinate to be ξ = ln(r + r0) leaves the mesh uniform in ξ, which is exactly what the geometry kernel already assumes -- an analytic function of the block corner and a constant spacing. The tree, the AMR level ladder, rnode, prolongation, coarsening and the ghost exchange all keep working unchanged in ξ; only the step from logical faces to physical ones takes an exp. It needs none of upstream MPI-AMRVAC's general stretched-grid machinery (qstretch, dxfirst, per-level stretch tables) -- which AGILE has removed -- and refinement and log stretching commute for free.

A positive log_r0 moves the singularity of the map from r = 0 out to r = −r0, making the mesh genuinely uniform (cell width r0·dξ) for r ≪ r0 and logarithmic for r ≫ r0. That is what lets a stretched radial grid reach r = 0: with the default log_r0 = 0 a log-radial grid can never reach the axis, so xprobmin1 ≤ 0 is rejected unless log_r0 > 0. logCylindrical with log_r0 > 0 therefore supports the cylindrical axis; logSpherical still has no usable origin (nothing fills the ghost cells across r = 0).

Two limitations of the stretching:

• The reconstruction stays uniform-stencil. The MUSCL reconstruction and the limiters assume equal cells, so the second-order constant degrades as the stretch ratio grows. Keep it modest -- the test cases run at ratios around 1.04 to 1.09 per cell. A positive log_r0 helps here, since the mesh is genuinely uniform in the region next to the axis where the stencil assumption would otherwise be worst.
• The timestep tightens. dr = (r + r0)·dξ, so a small log_r0 chosen for near-axis resolution buys a CFL restriction on top of the one the polar axis already imposes.

The polar axis

A spherical domain may span the full θ from 0 to π, and a cylindrical one may start at r = 0. The axis is not a boundary: a ghost cell just past the axis at azimuth φ is the interior cell at φ + π, so it is filled from the block half a revolution away (a "pi-periodic" treatment) rather than by a boundary condition. Setting it up takes three things in &meshlist/&boundlist, and check_pole_setup stops the run with a specific message if any is missing:

• 'pole' on the axis face(s) -- typeboundary_min2/max2 in spherical, typeboundary_min1 in cylindrical; it is all-or-nothing per face.
• φ periodic over exactly one full turn (xprobmin3 = 0, xprobmax3 = 1 in 2π units, periodic at both φ faces).
• an even number of level-1 blocks across φ.

The sign rule is the same in both systems: the component along the mirrored direction flips, the φ component flips, everything else is symmetric. Note that this makes the radial component odd at the cylindrical axis (m_r, B_r as well as m_φ, B_φ) -- upstream MPI-AMRVAC marks only the φ component antisymmetric there, which is wrong.

The polar axis is supported for hd, mhd, srhd and ffhd. Two things to know:

• The timestep collapses near the axis. ds(φ) = r·sin(θ)·dφ → 0 there, so setdt gives a severe CFL restriction. No axis filtering or ring averaging is implemented; keep pole cases at modest resolution.
• Volume averages cannot validate the pole copy -- a cell on the axis has vanishing volume and a zero-area axis face, so a wrong ghost value barely moves mean(w). The pole test cases therefore assert on the ghost cells directly against the analytic state at it = 0.

Current limitations

• AMR across curvilinear refinement levels is only lightly tested. Prolongation and coarsening volume-weight correctly, but fix_conserve is disabled for all geometries in src/mod_advance.fpp.
• General grid stretching is unavailable. A stretch_dim that varied the cell spacing within a block could not be expressed by the analytic metric derivation, so upstream MPI-AMRVAC's stretched-grid machinery was removed and the stretch_dim / qstretch_baselevel / nstretchedblocks_baselevel / stretch_uncentered keys are no longer read from &meshlist. The logarithmic radius above is the one supported non-uniform mesh, precisely because it is still uniform in ξ = ln(r + r0).
• divvector/curlvector in src/mod_geometry.fpp still assume a face lies midway between two cell centres, which is untrue for the volume barycentre. Neither routine is on a path this fork currently exercises.