Orbital Advection

Orbital Advection

Orbital advection is an algorithm to reduce the computational cost of hydrodynamic/MHD simulations in a shearing box or of astronomical disks (Masset 2000). This is also known as the FARGO (Fast Advection in Rotating Gaseous Objects) algorithm. By decomposing the velocity field v into the steady background orbital motion v_K and the residual v', the MHD equations have two parts: the typical MHD part and the linear orbital advection part. The orbital advection solves the equations separately using the operator splitting method. Since v_K is steady and known, the time step dt is determined regardless of the background orbital motion. As a result, large dt reduces the computational cost. Please see e.g., Stone & Gardinar 2010 for more details.

Configuration

Nothing to do at the configuration for orbital advection.

Input File

Parameter OAorder in the <orbital_advection> block is a flag for the orbital advection. By default, OAorder=0 is set and the orbital advection is turned off. Giving OAorder=1 activates the orbital advection using the first-order operator splitting method. Likewise, OAorder=2 specifies the second-order operator splitting method in the orbital advection. Using default orbital velocity profiles, there are other necessary parameters. In the Cartesian coordinates, those are qshear, the shearing rate (1.5 for Keplerian rotation), and local orbital frequency Omega0 as well as the shearing box. When the orbital advection and the shearing box are used together, shboxcoord = 1 is also required. Then, a typical <orbital_advection> block should look like:

    <orbital_advection>
    OAorder    = 2        # 0: w/o OA, 1: w/ OA(1st), 2: w/ OA(2nd)
    Omega0     = 1.0      # orbital frequency
    qshear     = 1.5      # shearing rate
    shboxcoord = 1        # 1: default x-y shear; 2: 2D x-z shear 

In the cylindrical and spherical polar coordinates, a positive value of the point source gravity GM in the <problem> block is needed. As an optional parameter, Omega0 in the <orbital_advection> block leads the rotating frame at the angular velocity Omega0 with the orbital advection. Then, the <orbital_advection> and <problem> blocks should be as follows:

    <orbital_advection>
    OAorder    = 2        # 0: w/o OA, 1: w/ OA(1st), 2: w/ OA(2nd)
    Omega0     = 1.0      # angular velocity of the rotating system (optional)
    <problem>
    GM         = 1.0      # point source gravity

In order to use the orbital advection, the orbital direction has to be periodic with uniform grid spacing. In Cartesian or cylindrical coordinates, the orbital direction is the x2 direction and one has to change the <mesh> block as follows:

    <mesh>
    ix2_bc     = periodic         # inner-X2 boundary flag
    ox2_bc     = periodic         # outer-X2 boundary flag
    x2rat      = 1.0              # X2 cell aspect ratio

In spherical polar coordinates, the orbital direction is the x3 direction and one has to change the <mesh> block as follows:

    <mesh>
    ix3_bc     = periodic         # inner-X3 boundary flag
    ox3_bc     = periodic         # outer-X3 boundary flag
    x3rat      = 1.0              # X3 cell aspect ratio

By default, the extra parameter orbital_system is set to false in the <output#> blocks, and the hydrodynamic variables in the outputs are based on v. To obtain the hydrodynamic variables based on v' rather than v, please set orbital_system = true.

    <output1>
    file_type      = vtk       # VTK data dump
    orbital_system = false     # false: default

Default Orbital Velocity Profile

Without a user-defined orbital velocity profile, the default orbital velocity profile is automatically set in OrbitalVelocityFunc OrbitalVelocity in the OrbitalAdvection class. Then, its derivatives in the directions perpendicular to the orbital direction are also set in OrbitalVelocityFunc OrbitalVelocityDerivative[i] (i = 0 or i = 1). One can access those using a pointer OrbitalAdvection *porb in the MeshBlock class. These functions are defined even if the orbital advection is not activated.

In Cartesian coordinates, v_K (x, y, z) = -q Ω₀ x. The value and the derivatives at x = (x, y, z) are obtained as follows:

  void MeshBlock::ProblemGenerator(ParameterInput *pin) {
    Real x, y, z;
    < set values in x, y, z >
    Real vel = porb->OrbitalVelocity(porb, x, y, z);              // vK(x, y, z)
    Real dv1 = porb->OrbitalVelocityDerivative[0](porb, x, y, z); // dvK/dx(x, y, z)
    Real dv2 = porb->OrbitalVelocityDerivative[1](porb, x, y, z); // dvK/dz(x, y, z)
  }

In cylindrical coordinates, v_K (R, φ, z) = sqrt(GM/R) - Ω₀ R for 2D and v_K (R, φ, z) = sqrt(GM/d) (R/d) - Ω₀ R for 3D, where d = sqrt(R² + z²). The value and the derivatives at x = (R, φ, z) are obtained as follows:

  void MeshBlock::ProblemGenerator(ParameterInput *pin) {
    Real R, phi, z;
    < set values in R, phi, z >
    Real vel = porb->OrbitalVelocity(porb, R, phi, z);              // vK(R, phi, z)
    Real dv1 = porb->OrbitalVelocityDerivative[0](porb, R, phi, z); // dvK/dR(R, phi, z)
    Real dv2 = porb->OrbitalVelocityDerivative[1](porb, R, phi, z); // dvK/dz(R, phi, z)
  }

In spherical polar coordinates, v_K (r, θ, φ) = sqrt(GM/(r sin(θ))) - Ω₀ r sin(θ). The value and the derivatives at x = (r, θ, φ) are obtained as follows:

  void MeshBlock::ProblemGenerator(ParameterInput *pin) {
    Real r, theta, phi;
    < set values in r, theta, phi >
    Real vel = porb->OrbitalVelocity(porb, r, theta, phi);              // vK(R, r, theta, phi)
    Real dv1 = porb->OrbitalVelocityDerivative[0](porb, r, theta, phi); // dvK/d(theta)(R, r, theta, phi)
    Real dv2 = porb->OrbitalVelocityDerivative[1](porb, r, theta, phi); // dvK/d(phi)(R, r, theta, phi)
  }

User-Defined Orbital Velocity Profile

In order to use a user-defined orbital velocity profile, they must be defined and then enrolled in Mesh::InitUserMeshData() in the problem generator.

First, the orbital velocity function must be defined and match the following function signature:

Real UserOrbitalVelocity(OrbitalAdvection * porb, Real x1, Real x2, Real x3);

The custom function can then be enrolled by calling:

EnrollOrbitalVelocity(UserOrbitalVelocity);

If using the user-defined functions for the derivatives of the orbital velocity, the functions must be defined an match the following function signature:

// x1 direction
Real UserOrbitalVelocityDerivative1(OrbitalAdvection * porb, Real x1, Real x2, Real x3);
// x3 direction in Cartesian and cylindrical, x2 direction in spherical polar
Real UserOrbitalVelocityDerivative2(OrbitalAdvection * porb, Real x1, Real x2, Real x3);

The custom functions can then be enrolled by calling:

// x1 direction
EnrollOrbitalVelocityDerivative(0, UserOrbitalVelocityDerivative1);
// x3 direction in Cartesian and cylindrical, x2 direction in spherical polar
EnrollOrbitalVelocityDerivative(1, UserOrbitalVelocityDerivative2);

Otherwise, the derivatives of the orbital velocity are calculated from the user-defined orbital velocity profile numerically.

Example

Example 1:

The 2D hydro shearing sheet test:

    > python configure.py --prob=ssheet --flux=hlle --eos=isothermal

and use input/hydro/athinput.ssheet.

Example 2:

The MHD shwave test:

    > python configure.py -b --prob=jgg --flux=hlld --eos=isothermal 

and use input/mhd/athinput.jgg

Example 3:

The 3D MRI test:

    > python configure.py -b --prob=hgb --flux=hlld --eos=isothermal 

and use input/mhd/athinput.hgb

Example 4:

The disk test in the cylindrical coordinates:

    > python configure.py --prob=disk --coord=cylindrical --flux=hlld

and use input/hydro/athinput.disk_cyl.

Example 5:

The disk test in the spherical polar coordinates:

    > python configure.py --prob=disk --coord=spherical_polar --flux=hlld

and use input/hydro/athinput.disk_sph.

Performance

Runs with the orbital advection require additional computational cost for each cycle in comparison with runs without the orbital advection. The cost comes from extra MPI communications and solving the linear advection. We find that the increase of the computational cost for one cycle is from 10% to 20% with the first-order splitting method and from 20% to 40% with the second-order splitting method. Therefore, it is expected that the orbital advection reduces the computational cost when the orbital speed is faster by 20% for OAorder=1 and by 40% OAorder=2 than the sound speed.

Note

• There are two restrictions on the time step dt from the orbital advection. One requires that the shear distance between cells used in the reconstruction at each cell is smaller than the half cell size in the orbital direction for the first-order scheme and than the cell size in the orbital direction for the second-order scheme. The other prohibits the orbital motion for one cycle from exceeding the meshblock size in the orbital direction.
• The orbital velocity profile must not change in time.
• Orbital advection works with mesh refinement (AMR and SMR), passive scalars, shearing box, diffusion processes (hydro and field), general EOS, self-gravity (FFT and multigrid), and super time stepping.
• The orbital advection works with neither special relativity nor general relativity.
• Both the PLM reconstructions (time/xorder=2 or time/xorder=2c) and the PPM reconstructions (time/xorder=3 or time/xorder=3c) are available with the orbital advection. When giving time/xorder=4, however, the orbital remapping will be performed with only the 3rd-order spatial accuracy like the shear-periodic boundaries.
• Neither time/integrator=rk4 nor time/integrator=ssprk5_4 are allowed with both the orbital advection and the mesh refinement.
• When using orbital advection with STS, orbital_advection/OAorder=2 and time/sts_integrator=rkl2 are recommended.