PrincetonUniversity/athena

Implementing `gr_torus` for user defined spacetime

osokoliuk opened this issue · 1 comments

Hi. I'm quite new to Athena++, so sorry if this question is confusing . I'm trying to somehow adapt GR MHD torus for user defined spacetime, admitting spherical symmetry. For example, I have a spacetime (Schild wormhole with throat radius r0 and present tidal forces):

#include <cmath>  // abs(), cos(), log(), sin(), sqrt()


#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "../eos/eos.hpp"
#include "../mesh/mesh.hpp"
#include "../parameter_input.hpp"
#include "coordinates.hpp"

void SchildWH(Real x1, Real x2, Real x3, ParameterInput *pin,
    AthenaArray<Real> &g, AthenaArray<Real> &g_inv, AthenaArray<Real> &dg_dx1,
    AthenaArray<Real> &dg_dx2, AthenaArray<Real> &dg_dx3)
{
  // Extract inputs
  Real aa = pin->GetReal("coord", "aa");
  Real r0 = pin->GetReal("coord", "r0");

  Real r = x1;
  Real theta = x2;
  Real phi = x3;

  // Calculate intermediate quantities
  Real sth = std::sin(theta);
  Real cth = std::cos(theta);
  Real shape = r0;
  Real redshift = -aa/r;
  Real dshape = 0;
  Real dredshift =  aa/(r*r);
  
  // Set covariant components
  g(I00) = -std::exp(2*redshift);
  g(I11) = 1/(1-shape/r);
  g(I22) = r*r;
  g(I33) = r*r*sth*sth;
  g(I01) = 0.0;
  g(I02) = 0.0;
  g(I03) = 0.0;
  g(I12) = 0.0;
  g(I13) = 0.0;
  g(I23) = 0.0;

  // Set contravariant components
  g_inv(I00) = -std::exp(-2*redshift);
  g_inv(I11) = 1-shape/r;
  g_inv(I22) = 1/(r*r);
  g_inv(I33) = 1/(r*r*sth*sth);
  g_inv(I01) = 0.0;
  g_inv(I02) = 0.0;
  g_inv(I03) = 0.0;
  g_inv(I12) = 0.0;
  g_inv(I13) = 0.0;
  g_inv(I23) = 0.0;

  // Set r-derivatives of covariant components
  dg_dx1(I00) = -2*std::exp(-2*redshift)*dredshift;
  dg_dx1(I11) = (r*dshape-shape)/((r-shape)*(r-shape));
  dg_dx1(I22) = 2*r;
  dg_dx1(I33) = 2*r*sth*sth
  dg_dx1(I01) = 0.0;
  dg_dx1(I02) = 0.0;
  dg_dx1(I03) = 0.0;
  dg_dx1(I12) = 0.0;
  dg_dx1(I13) = 0.0;
  dg_dx1(I23) = 0.0;

  // Set theta-derivatives of covariant components
  dg_dx2(I00) = 0;
  dg_dx2(I11) = 0;
  dg_dx2(I22) = 0;
  dg_dx2(I33) = 2*r*r*sth*cth
  dg_dx2(I01) = 0.0;
  dg_dx2(I02) = 0.0;
  dg_dx2(I03) = 0.0;
  dg_dx2(I12) = 0.0;
  dg_dx2(I13) = 0.0;
  dg_dx2(I23) = 0.0;

  // Set phi-derivatives of covariant components
  for (int n = 0; n < NMETRIC; ++n) {
    dg_dx3(n) = 0.0;
  }
  return;
}

Is it possible to adapt this metric (or any other gr_user metric) for GR torus, or it works only for KerrSchild case?

The gr_torus.cpp file very much assumes a Kerr spacetime. You might be able to augment the coordinate transformation functions like GetKerrSchildCoordinates and TransformContravariantFromBoyerLindquist with a new case to get it to produce something, but a more fundamental issue is whether this would be meaningful. The file implements a hydrostatic solution particular to Kerr, and a different spacetime will fail to achieve this property if the same prescription is followed.

It might be worth writing a new (and much shorter) problem generator that initializes the fluid in a simple, near-equilibrium fashion.