Launching jets via a source term in SR

[paytonrodman]

Hi all!

I am trying to simulate jets in a cartesian grid and so am using the EnrollUserExplicitSourceFunction to set relevant jet variables (e.g. density, pressure, velocities) within a couple of cylindrical regions inside the mesh, rather than at a boundary (example pictures below).

I can set these values fine, but the simulation isn't evolving - there's a non-zero vy_jet inside the launching cylinders but nothing leaves that region and all variables are static.

The code is based heavily on the jet.cpp example with SR calculations from the gr_shock_tube.cpp example.

I'm not sure if it's a conflict between SR and the source term, or if perhaps the sharp discontinuities on the launch regions are messing with energy conservation.

Any help is greatly appreciated!

Your Athena++ distribution has now been configured with the following options:
  Problem generator:          jet
  Coordinate system:          cartesian
  Equation of state:          adiabatic
  Riemann solver:             hllc
  Magnetic fields:            OFF
  Number of scalars:          0
  Special relativity:         ON
  General relativity:         OFF
  Frame transformations:      OFF
  Self-Gravity:               OFF
  Super-Time-Stepping:        OFF
  Debug flags:                OFF
  Code coverage flags:        OFF
  Linker flags:                -L/usr/local/Cluster-Apps/hdf5/mpi/intel/2019.3/1.10.5/lib  -lhdf5
  Floating-point precision:   double
  Number of ghost cells:      2
  MPI parallelism:            ON
  OpenMP parallelism:         OFF
  FFT:                        OFF
  HDF5 output:                ON
  HDF5 precision:             single
  Compiler:                   icpc

//----------------------------------------------------------------------------------------
//! \f  Total source term

void Total_SrcTerm(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) {

  //===================================
  // Add jets
  //===================================
  // Prepare auxiliary array
  AthenaArray<Real> bb_src;
  bb_src.NewAthenaArray(3, pmb->ke+1, pmb->je+1, pmb->ie+1);

  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) {
        Real vy_jet_local;
        Real ycoord = pmb->pcoord->x2v(j);
        Real rad = std::sqrt(SQR(pmb->pcoord->x1v(i)) + SQR(pmb->pcoord->x3v(k)));

        // Check if within jet launching cylinders
        if (rad <= r_jet) {
          if (std::abs(ycoord) > ymin_jet && std::abs(ycoord) < ymax_jet) {

            // Set negative values for lower jet
            if (ycoord < 0.) {
              vy_jet_local = -vy_jet;
            } else {
              vy_jet_local = vy_jet;
            }

            // Construct 4-vectors
            Real ut = std::sqrt(1.0 / (1.0 - (SQR(vx_jet)+SQR(vy_jet_local)+SQR(vz_jet))));
            Real ux = ut * vx_jet;
            Real uy = ut * vy_jet_local;
            Real uz = ut * vz_jet;
            Real bt = bx_jet*ux + by_jet*uy + bz_jet*uz;
            Real bx = (bx_jet + bt * ux) / ut;
            Real by = (by_jet + bt * uy) / ut;
            Real bz = (bz_jet + bt * uz) / ut;

            // Set hydro primitives
            pmb->phydro->w(IDN,k,j,i) = pmb->phydro->w1(IDN,k,j,i) = d_jet;
            pmb->phydro->w(IPR,k,j,i) = pmb->phydro->w1(IPR,k,j,i) = p_jet;
            pmb->phydro->w(IVX,k,j,i) = pmb->phydro->w1(IM1,k,j,i) = ux;
            pmb->phydro->w(IVY,k,j,i) = pmb->phydro->w1(IM2,k,j,i) = uy;
            pmb->phydro->w(IVZ,k,j,i) = pmb->phydro->w1(IM3,k,j,i) = uz;

            // Define cell-centered magnetic fields
            bb_src(IB1,k,j,i) = bx * ut - bt * ux;
            bb_src(IB2,k,j,i) = by * ut - bt * uy;
            bb_src(IB3,k,j,i) = bz * ut - bt * uz;
          }
        }
      }
    }
  }
  // Update conserved variables
  pmb->peos->PrimitiveToConserved(pmb->phydro->w, bb_src, pmb->phydro->u, pmb->pcoord,
                                  pmb->is, pmb->ie, pmb->js, pmb->je, pmb->ks, pmb->ke);

  // set magnetic field in jet inlet zones
  if (MAGNETIC_FIELDS_ENABLED) {
    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) {
          Real ycoord = pmb->pcoord->x2v(j);
          Real rad = std::sqrt(SQR(pmb->pcoord->x1v(i)) + SQR(pmb->pcoord->x3v(k)));

          if (rad <= r_jet) {
            if (std::abs(ycoord) > ymin_jet && std::abs(ycoord) < ymax_jet) {
              pmb->pfield->b.x1f(k,j,i) = bx_jet;
              pmb->pfield->b.x2f(k,j,i) = by_jet;
              pmb->pfield->b.x3f(k,j,i) = bz_jet;
            }
          }
        }
      }
    }
  }
}