Skip to content

Commit

Permalink
Merge branch 'main' into c2p
Browse files Browse the repository at this point in the history
  • Loading branch information
@jaykalinani
jaykalinani committed Aug 13, 2024
2 parents eebc627 + 085dd8a commit bc96a82
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 44 deletions.
4 changes: 4 additions & 0 deletions AsterX/param.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,10 @@ USES CCTK_REAL rho_strict
USES BOOLEAN Ye_lenient
USES CCTK_INT max_iter
USES CCTK_REAL c2p_tol
USES CCTK_REAL r_atmo
USES CCTK_REAL n_rho_atmo
USES CCTK_REAL n_press_atmo
USES BOOLEAN thermal_eos_atmo

SHARES: ReconX
USES KEYWORD reconstruction_method
Expand Down
51 changes: 34 additions & 17 deletions AsterX/src/con2prim.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,28 +55,45 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,

const smat<GF3D2<const CCTK_REAL>, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};

// Setting up atmosphere
const CCTK_REAL rho_atmo_cut = rho_abs_min * (1 + atmo_tol);
const CCTK_REAL gm1 = eos_cold.gm1_from_valid_rmd(rho_abs_min);
CCTK_REAL eps_atm = eos_cold.sed_from_valid_gm1(gm1);
eps_atm = std::min(std::max(eos_th.rgeps.min, eps_atm), eos_th.rgeps.max);
const CCTK_REAL p_atm =
eos_th.press_from_valid_rho_eps_ye(rho_abs_min, eps_atm, Ye_atmo);
atmosphere atmo(rho_abs_min, eps_atm, Ye_atmo, p_atm, rho_atmo_cut);

// Construct Noble c2p object:
c2p_2DNoble c2p_Noble(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
B_lim, Ye_lenient);

// Construct Palenzuela c2p object:
c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
B_lim, Ye_lenient);

// Loop over the interior of the grid
cctk_grid.loop_int_device<
1, 1, 1>(grid.nghostzones, [=] CCTK_DEVICE(
const PointDesc
&p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
// Setting up atmosphere
CCTK_REAL rho_atm = 0.0; // dummy initialization
CCTK_REAL press_atm = 0.0; // dummy initialization
CCTK_REAL eps_atm = 0.0; // dummy initialization
CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);

// Grading rho
rho_atm = (radial_distance > r_atmo)
? (rho_abs_min * pow((r_atmo / radial_distance), n_rho_atmo))
: rho_abs_min;
const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol);

// Grading pressure based on either cold or thermal EOS
if (thermal_eos_atmo) {
press_atm = (radial_distance > r_atmo)
? (p_atmo * pow(r_atmo / radial_distance, n_press_atmo))
: p_atmo;
eps_atm = eos_th.eps_from_valid_rho_press_ye(rho_atm, press_atm, Ye_atmo);
} else {
const CCTK_REAL gm1 = eos_cold.gm1_from_valid_rmd(rho_atm);
eps_atm = eos_cold.sed_from_valid_gm1(gm1);
eps_atm = std::min(std::max(eos_th.rgeps.min, eps_atm), eos_th.rgeps.max);
press_atm = eos_th.press_from_valid_rho_eps_ye(rho_atm, eps_atm, Ye_atmo);
}
atmosphere atmo(rho_atm, eps_atm, Ye_atmo, press_atm, rho_atmo_cut);

// Construct Noble c2p object:
c2p_2DNoble c2p_Noble(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
B_lim, Ye_lenient);

// Construct Palenzuela c2p object:
c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, rho_strict,
vw_lim, B_lim, Ye_lenient);

/* Get covariant metric */
const smat<CCTK_REAL, 3> glo(
[&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); });
Expand Down
74 changes: 47 additions & 27 deletions AsterX/src/con2prim_rpa.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,37 +25,57 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {

const smat<GF3D2<const CCTK_REAL>, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};

// Setting up initial data EOS
const CCTK_REAL n = 1 / (poly_gamma - 1); // Polytropic index
const CCTK_REAL adiab_ind_id = 1.0 / (poly_gamma - 1);
const CCTK_REAL rmd_p = pow(poly_k, -n); //Polytropic density scale
const auto eos_id = make_eos_barotr_poly(adiab_ind_id, rmd_p, rho_max);

// Setting up evolution EOS
const CCTK_REAL adiab_ind_evol = 1.0 / (gl_gamma - 1);
const auto eos = make_eos_idealgas(adiab_ind_evol, eps_max, rho_max);

// Setting up atmosphere
const CCTK_REAL rho_atmo_cut = rho_abs_min * (1 + atmo_tol);
CCTK_REAL eps_atm = eos_id.at_rho(rho_abs_min).eps();
eps_atm = eos.range_eps(rho_abs_min, Ye_atmo).limit_to(eps_atm);
CCTK_REAL p_atm = eos.at_rho_eps_ye(rho_abs_min, eps_atm, Ye_atmo).press();
const atmosphere atmo(rho_abs_min, eps_atm, Ye_atmo, p_atm, rho_atmo_cut);

CCTK_REAL dummy_Ye = 0.5;
CCTK_REAL dummy_dYe = 0.5;

// Get a recovery function
con2prim_mhd cv2pv(eos, 1e-5, 1, 100, 100, atmo, 1e-8, 100);
// con2prim_mhd cv2pv(eos, rho_strict, ye_lenient, max_z, max_b, atmo,
// c2p_acc,
// max_iter);

// Loop over the interior of the grid
cctk_grid.loop_int_device<
1, 1, 1>(grid.nghostzones, [=] CCTK_DEVICE(
const PointDesc
&p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
// Setting up initial data EOS
const CCTK_REAL n = 1 / (poly_gamma - 1); // Polytropic index
const CCTK_REAL adiab_ind_id = 1.0 / (poly_gamma - 1);
const CCTK_REAL rmd_p = pow(poly_k, -n); // Polytropic density scale
const auto eos_id = make_eos_barotr_poly(adiab_ind_id, rmd_p, rho_max);

// Setting up evolution EOS
const CCTK_REAL adiab_ind_evol = 1.0 / (gl_gamma - 1);
const auto eos = make_eos_idealgas(adiab_ind_evol, eps_max, rho_max);

// Setting up atmosphere

CCTK_REAL rho_atm = 0.0; // dummy initialization
CCTK_REAL press_atm = 0.0; // dummy initialization
CCTK_REAL eps_atm = 0.0; // dummy initialization
CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);

// Grading rho
rho_atm = (radial_distance > r_atmo)
? (rho_abs_min * pow((r_atmo / radial_distance), n_rho_atmo))
: rho_abs_min;
const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol);

// Grading pressure and eps
if (thermal_eos_atmo) {
press_atm = (radial_distance > r_atmo)
? (p_atmo * pow(r_atmo / radial_distance, n_press_atmo))
: p_atmo;
//TODO: eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps() does not exist in RePrimAnd
//Currently computing eps from ideal gas EOS
//eps_atm = eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps();
eps_atm = press_atm/(rho_atm*(gl_gamma - 1.));
} else {
eps_atm = eos_id.at_rho(rho_atm).eps();
eps_atm = eos.range_eps(rho_atm, Ye_atmo).limit_to(eps_atm);
press_atm = eos.at_rho_eps_ye(rho_atm, eps_atm, Ye_atmo).press();
}
const atmosphere atmo(rho_atm, eps_atm, Ye_atmo, press_atm, rho_atmo_cut);

CCTK_REAL dummy_Ye = 0.5;
CCTK_REAL dummy_dYe = 0.5;

// Get a recovery function
con2prim_mhd cv2pv(eos, rho_strict, Ye_lenient, vw_lim, B_lim, atmo, c2p_tol,
max_iter);

/* Get covariant metric */
const smat<CCTK_REAL, 3> glo(
[&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); });
Expand Down Expand Up @@ -117,7 +137,7 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
// velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I),
Avec_x(p.I), Avec_y(p.I), Avec_z(p.I));
}
//set to atmo
// set to atmo
cv.bcons(0) = dBx(p.I);
cv.bcons(1) = dBy(p.I);
cv.bcons(2) = dBz(p.I);
Expand Down
22 changes: 22 additions & 0 deletions Con2PrimFactory/param.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,3 +70,25 @@ CCTK_REAL c2p_tol "c2p torelance for root finding" STEERABLE=ALWAYS
{
0:* :: "Larger than zero"
} 1e-10

# Parameters for graded atmosphere

CCTK_REAL r_atmo "Radial distance up to which rho_abs_min is kept constant, and after which it decreases as a power law function." STEERABLE=ALWAYS
{
(0.0:* :: ""
} 1.0e100

CCTK_REAL n_rho_atmo "Exponential of the radial power law atmosphere prescription for density." STEERABLE=ALWAYS
{
0.0:* :: ""
} 0.0

CCTK_REAL n_press_atmo "Exponential of the radial power law atmosphere prescription for pressure." STEERABLE=ALWAYS
{
0.0:* :: ""
} 0.0

BOOLEAN thermal_eos_atmo "Whether to use the thermal (evolution) EOS for setting pressure in the atmosphere. If set to no, cold (initial data) EOS is used instead." STEERABLE=ALWAYS
{
} "no"

0 comments on commit bc96a82

Please sign in to comment.