Skip to content

[WIP] Radiative GRMHD for KHARMA #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: dev
Choose a base branch
from
Open

[WIP] Radiative GRMHD for KHARMA #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,4 +33,7 @@ include_directories(external/parthenon/src)
include_directories(external/variant/include)

# KHARMA folder
add_subdirectory(kharma)
#add_subdirectory(kharma)

# KIPOLE folder
add_subdirectory(kipole)
149 changes: 79 additions & 70 deletions kharma/coordinates/coordinate_systems.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,69 @@ using mpark = std;
using namespace parthenon;
using GReal = Real;

template<typename Function>
KOKKOS_FUNCTION void root_find(const GReal Xembed[GR_DIM], GReal Xnative[GR_DIM], Function coord_to_embed);
/**
* Root finder for X[2] since it is sometimes not analytically invertible
* Written so it can be extended if people have very crazy coordinate ideas that require 2D solves
* @param Xembed the vector of embedding coordinates to convert
* @param Xnative output; but should have all native coordinates except X[2] already
* @param coord_to_embed function taking the vector X to embedding coordinates
*/
template<typename T>
KOKKOS_INLINE_FUNCTION void root_find(const GReal Xembed[GR_DIM], GReal Xnative[GR_DIM], T transform)
{
double th = Xembed[2];
double tha, thb, thc;

// Currently only solves in X[2] but could be multi-dimensional
double Xa[GR_DIM], Xb[GR_DIM], Xc[GR_DIM], Xtmp[GR_DIM];
Xa[1] = Xnative[1];
Xa[3] = Xnative[3];

Xb[1] = Xa[1];
Xb[3] = Xa[3];
Xc[1] = Xa[1];
Xc[3] = Xa[3];

if (Xembed[2] < M_PI / 2.) {
Xa[2] = 0.;
Xb[2] = 0.5 + SINGSMALL;
} else {
Xa[2] = 0.5 - SINGSMALL;
Xb[2] = 1.;
}

double tol = 1.e-9;
transform->coord_to_embed(Xa, Xtmp);
tha = Xtmp[2];
transform->coord_to_embed(Xb, Xtmp);
thb = Xtmp[2];

// check limits first
if (fabs(tha-th) < tol) {
Xnative[2] = Xa[2];
return;
} else if (fabs(thb-th) < tol) {
Xnative[2] = Xb[2];
return;
}

// bisect for a bit
for (int i = 0; i < 1000; i++) {
Xc[2] = 0.5 * (Xa[2] + Xb[2]);
transform->coord_to_embed(Xc, Xtmp);
thc = Xtmp[2];

if ((thc - th) * (thb - th) < 0.)
Xa[2] = Xc[2];
else
Xb[2] = Xc[2];

if (fabs(thc - th) < tol)
break;
}

Xnative[2] = Xc[2];
}

/**
* EMBEDDING SYSTEMS:
Expand Down Expand Up @@ -123,11 +184,17 @@ class SphKSCoords {
DLOOP2 vcon[mu] += trans[mu][nu]*vcon_bl[nu];
}

// TODO more: isco etc.
// TODO more functions?
KOKKOS_INLINE_FUNCTION GReal rhor() const
{
return (1. + sqrt(1. - a*a));
}
KOKKOS_INLINE_FUNCTION GReal r_isco() const
{
double z1 = 1. + pow(1. - a * a, 1. / 3.) * (pow(1. + a, 1. / 3.) + pow(1. - a, 1. / 3.));
double z2 = sqrt(3. * a * a + z1 * z1);
return 3. + z2 - copysign(sqrt((3. - z1) * (3. + z1 + 2. * z2)), a);
}
};

/**
Expand Down Expand Up @@ -175,6 +242,12 @@ class SphBLCoords {
{
return (1. + sqrt(1. - a*a));
}
KOKKOS_INLINE_FUNCTION GReal r_isco() const
{
double z1 = 1. + pow(1. - a * a, 1. / 3.) * (pow(1. + a, 1. / 3.) + pow(1. - a, 1. / 3.));
double z2 = sqrt(3. * a * a + z1 * z1);
return 3. + z2 - copysign(sqrt((3. - z1) * (3. + z1 + 2. * z2)), a);
}
};

/**
Expand Down Expand Up @@ -286,7 +359,7 @@ class ModifyTransform {
Xnative[1] = log(Xembed[1]);
Xnative[3] = Xembed[3];
// Treat the special case
root_find(Xembed, Xnative, &ModifyTransform::coord_to_embed);
root_find(Xembed, Xnative, this);
}
/**
* Transformation matrix for contravariant vectors to embedding, or covariant vectors to native
Expand Down Expand Up @@ -363,7 +436,7 @@ class FunkyTransform {
Xnative[1] = log(Xembed[1]);
Xnative[3] = Xembed[3];
// Treat the special case
root_find(Xembed, Xnative, &FunkyTransform::coord_to_embed);
root_find(Xembed, Xnative, this);
}
/**
* Transformation matrix for contravariant vectors to embedding, or covariant vectors to native
Expand Down Expand Up @@ -410,68 +483,4 @@ class FunkyTransform {
// Bundle coordinates and transforms into umbrella variant types
// Note nesting isn't allowed -- do it yourself by calling the steps if that's really important...
using SomeBaseCoords = mpark::variant<SphMinkowskiCoords, CartMinkowskiCoords, SphBLCoords, SphKSCoords>;
using SomeTransform = mpark::variant<SphNullTransform, CartNullTransform, ModifyTransform, FunkyTransform>;

/**
* Root finder for X[2] since it is sometimes not analytically invertible
* Written so it can be extended if people have very crazy coordinate ideas that require 2D solves
* @param Xembed the vector of embedding coordinates to convert
* @param Xnative output; but should have all native coordinates except X[2] already
* @param coord_to_embed function taking the vector X to embedding coordinates
*/
template<typename Function>
KOKKOS_FUNCTION void root_find(const GReal Xembed[GR_DIM], GReal Xnative[GR_DIM], Function& coord_to_embed)
{
double th = Xembed[2];
double tha, thb, thc;

// Currently only solves in X[2] but could be multi-dimensional
double Xa[GR_DIM], Xb[GR_DIM], Xc[GR_DIM], Xtmp[GR_DIM];
Xa[1] = Xnative[1];
Xa[3] = Xnative[3];

Xb[1] = Xa[1];
Xb[3] = Xa[3];
Xc[1] = Xa[1];
Xc[3] = Xa[3];

if (Xembed[2] < M_PI / 2.) {
Xa[2] = 0.;
Xb[2] = 0.5 + SINGSMALL;
} else {
Xa[2] = 0.5 - SINGSMALL;
Xb[2] = 1.;
}

double tol = 1.e-9;
coord_to_embed(Xa, Xtmp);
tha = Xtmp[2];
coord_to_embed(Xb, Xtmp);
thb = Xtmp[2];

// check limits first
if (fabs(tha-th) < tol) {
Xnative[2] = Xa[2];
return;
} else if (fabs(thb-th) < tol) {
Xnative[2] = Xb[2];
return;
}

// bisect for a bit
for (int i = 0; i < 1000; i++) {
Xc[2] = 0.5 * (Xa[2] + Xb[2]);
coord_to_embed(Xc, Xtmp);
thc = Xtmp[2];

if ((thc - th) * (thb - th) < 0.)
Xa[2] = Xc[2];
else
Xb[2] = Xc[2];

if (fabs(thc - th) < tol)
break;
}

Xnative[2] = Xc[2];
}
using SomeTransform = mpark::variant<SphNullTransform, CartNullTransform, ModifyTransform, FunkyTransform>;
Loading