Adding constant momentum source in cylindrical coordinates

[rahulgaur104]

Hi everyone,

I want to add a constant force/volume(force density?) source in the azimuthal direction to the cylindrical system(r, phi, z) that I am simulating. The snippet of my code that does this is shown below:

//----------------------------------------------------------------------------------------
//! User-defined azimuthal momentum source term throughout the plasma
//----------------------------------------------------------------------------------------

void ExtSource(MeshBlock *pmb, const Real time, const Real dt,
              const AthenaArray<Real> &prim, const AthenaArray<Real> &prim_scalar,
              const AthenaArray<Real> &bcc, AthenaArray<Real> &cons,
              AthenaArray<Real> &cons_scalar) {


  for (int k=pmb->ks; k<=pmb->ke; ++k) {
    for (int j=pmb->js; j<=pmb->je; ++j) {
      for (int i=pmb->is; i<=pmb->ie; ++i) {
         cons(IM2,k,j,i) += SourceStrength*110000*dt*(1/(1 + std::exp(-60*(time-0.01)))) ;
        }
      }
    }
  return;
}

where SourceStrength is a number that is provided in the athinput file. Is this the right way to add a force/volume source in cylindrical coordinates that is independent of x1? I was looking at Skinner and Ostriker's paper (here) and I saw the source term in the momentum equation that contained factors of x1(or r).

To understand this, I looked at line 320 of /src/coordinates/cylindrical.cpp which confused me even further. Is coord_src2_i_(i) the same as the constant momentum density source term SourceStrength?

My ultimate objective is to add a constant force/volume source in the azimuthal(x2(or phi)) direction. Any feedback would be appreciated.

Thank you,
Rahul

[c-white]

The code you wrote should be correct.

The paper you cite (see also the equations near 2.14 of this PhD thesis, for more details) is adding a phi-momentum source term, but in a different context. The conservation laws/evolution equations have the form

∂_t(conserved) + div(flux) = S_physical(primitive).

In non-Cartesian coordinates, the divergence operator has extra/modified terms, so this becomes

∂_t(conserved) + ∂_x1(flux1) + ∂_x2(flux2) + ∂_x3(flux3) = S_geometric(primitive) + S_physical(primitive).

The problem is, if you naively integrate this equation over a single cell to turn it into a finite-volume scheme, the fluxes evaluated at cell faces won't combine with the cell-averaged geometric source term to exactly get the discrete divergence operator you want, causing things like secular errors in angular momentum. As it turns out, you can write the geometric source term in terms of face-centered fluxes (rather than cell-centered conserved or primitive variables):

∂_t(conserved) + ∂_x1(flux1) + ∂_x2(flux2) + ∂_x3(flux3) = S'_geometric(fluxes) + S_physical(primitive).

Even if S_geometric doesn't depend on position, S'_geometric may well do so.

However, none of this affects the physical source terms you are likely to be implementing as a user.

[rahulgaur104]

Thanks, Chris! If I understand you correctly, the user-implemented term SourceStrength*dt, SourceStrength being a constant, has dimensions of Force/volume and is independent of all coordinates.

[c-white]

Yes, SourceStrength should be force/volume.