Implementing gr_torus for user defined spacetime

[osokoliuk]

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?

[c-white]

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.

[osokoliuk]

I'll try to write a new problem generator for a near-equilibrium fluid accretion disk. However, I'm not really a skilled programmer, so it will be quite a challenge. Thank you for your fast reply!

[osokoliuk]

I managed to adapt Bondi spherical accretion to the wormhole spacetime. But, when I'm trying to plot the results with

~/athena/vis/python/plot_spherical.py SchildWH.out1.00001.athdf rho show

I get an error:

Traceback (most recent call last):
  File "/home/oleksii/athena/vis/python/plot_spherical.py", line 395, in <module>
    main(**vars(args))
  File "/home/oleksii/athena/vis/python/plot_spherical.py", line 82, in main
    data = athena_read.athdf(kwargs['data_file'], quantities=quantities,
  File "/home/oleksii/athena/vis/python/athena_read.py", line 505, in athdf
    raise AthenaError('Coordinates not recognized')
athena_read.AthenaError: Coordinates not recognized

How could this problem be solved? I'm attaching the pgen file and input file:
https://www.mediafire.com/file/3rkjilpxa3jkg5s/SchildWH.cpp/file
https://www.mediafire.com/file/4fcmvj730ptrv6l/athinput.SchildWH_inputs/file

Configure them with

./configure.py -g --prob SchildWH --coord=gr_user --flux hlle --nghost 4 -hdf5

[c-white]

The plotting scripts call the athdf function in vis/python/athena_read.py. That function is concerned with assembling the data into a numpy array as accurately as possible, doing the right thing even in cases with mesh refinement or non-uniform spacing. Sometimes this involves inferring coordinate data not explicitly saved in the data dump, and for that it needs to know which of the many coordinate systems were being used. Unfortunately, with user-specified coordinates it doesn't really know what to do.

In this particular case, it can't figure out how to calculate the cell-centered coordinates given the coordinates of the faces. If you were calling athdf directly, you could provide the optional center_func_* functions as inputs (see the documentation). Currently, however, there's no way to pass this through the plotting script.

Instead, I recommend modifying athena_read.py in 3 places. Just add in or coord == 'gr_user' to lines 492, 509, and 528. This way, the code will assume the cell center is halfway between the two faces. Actually this is probably always the case for the code proper, so I should probably make this change part of the code base.

[osokoliuk]

Thank you, that solved an issue. I have the last question - how could I point out exactly what gr_user spacetime I want to use? Even if I defined it and enrolled in the Mesh::InitUserMeshData, Bondi accretion pgen that I wrote still needs the values of mass for Schwarzschild black hole in the inputs

[c-white]

If your pgen needs a parameter, you can declare it near the the top of the file, read in its value in InitUserMeshData, and then use the variable in ProblemGenerator or any auxiliary function you may have defined. See for example gamma_adi in the existing gr_torus.cpp file.