Skip to content

Commit

Permalink
AsterX: Add code in C2P for BH interior handling
Browse files Browse the repository at this point in the history
  • Loading branch information
@jaykalinani
jaykalinani committed Aug 27, 2024
1 parent 43fdc63 commit 9e28ade
Show file tree
Hide file tree
Showing 4 changed files with 194 additions and 12 deletions.
4 changes: 4 additions & 0 deletions AsterX/param.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,10 @@ USES CCTK_REAL r_atmo
USES CCTK_REAL n_rho_atmo
USES CCTK_REAL n_press_atmo
USES BOOLEAN thermal_eos_atmo
USES CCTK_REAL alp_thresh
USES CCTK_REAL rho_BH
USES CCTK_REAL eps_BH
USES CCTK_REAL vwlim_BH

SHARES: ReconX
USES KEYWORD reconstruction_method
Expand Down
1 change: 1 addition & 0 deletions AsterX/schedule.ccl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,7 @@ SCHEDULE AsterX_Con2Prim IN AsterX_Con2PrimGroup AFTER AsterX_ComputedBFromdBsta
{
LANG: C
READS: ADMBaseX::metric(interior)
READS: ADMBaseX::lapse(everywhere)
READS: dens(interior) tau(interior) mom(interior) dB(interior)
READS: saved_prims(interior)
READS: Avec_x(interior) Avec_y(interior) Avec_z(interior)
Expand Down
67 changes: 63 additions & 4 deletions AsterX/src/con2prim.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,33 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,
atmo.set(pv_seeds);
}

// Modifying primitive seeds within BH interiors before C2Ps are called
// NOTE: By default, alp_thresh=0 so the if condition below is never
// triggered. One must be very careful when using this functionality and
// must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the parfile

if (alp(p.I) < alp_thresh) {
if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) {
pv_seeds.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of
// initial NS or disk
pv_seeds.eps = eps_BH;
pv_seeds.Ye = Ye_atmo;
pv_seeds.press =
eos_th.press_from_valid_rho_eps_ye(rho_BH, eps_BH, Ye_atmo);
if ( ((pv_seeds.vel(0) * pv_seeds.w_lor) > vwlim_BH) ||
((pv_seeds.vel(1)*pv_seeds.w_lor) > vwlim_BH) ||
((pv_seeds.vel(2)*pv_seeds.w_lor) > vwlim_BH) ) {
CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH);
CCTK_REAL vlim_BH = vwlim_BH / wlim_BH;
pv_seeds.vel(0) *= vlim_BH / pv_seeds.vel(0);
pv_seeds.vel(1) *= vlim_BH / pv_seeds.vel(1);
pv_seeds.vel(2) *= vlim_BH / pv_seeds.vel(2);
pv_seeds.w_lor = wlim_BH;
}
cv.from_prim(pv_seeds, glo);
}
}

// Construct error report object:
c2p_report rep_first;
c2p_report rep_second;
Expand Down Expand Up @@ -169,6 +196,35 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,
if (rep_first.failed() && rep_second.failed()) {
printf("Second C2P failed too :( :( \n");
rep_second.debug_message();

// Treatment for BH interiors after C2P failures
// NOTE: By default, alp_thresh=0 so the if condition below is never
// triggered. One must be very careful when using this functionality and
// must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the
// parfile
if (alp(p.I) < alp_thresh) {
if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) {
pv.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of initial
// NS or disk
pv.eps = eps_BH;
pv.Ye = Ye_atmo;
pv.press =
eos_th.press_from_valid_rho_eps_ye(rho_BH, eps_BH, Ye_atmo);
if (((pv.vel(0) * pv.w_lor) > vwlim_BH) ||
((pv.vel(1)*pv.w_lor) > vwlim_BH) ||
((pv.vel(2)*pv.w_lor) > vwlim_BH) ) {
CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH);
CCTK_REAL vlim_BH = vwlim_BH / wlim_BH;
pv.vel(0) *= vlim_BH / pv.vel(0);
pv.vel(1) *= vlim_BH / pv.vel(1);
pv.vel(2) *= vlim_BH / pv.vel(2);
pv.w_lor = wlim_BH;
}
cv.from_prim(pv, glo);
rep_first.set_atmo = 0;
rep_second.set_atmo = 0;
}
}
con2prim_flag(p.I) = 0;
}

Expand Down Expand Up @@ -218,13 +274,16 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,
pv.scatter(rho(p.I), eps(p.I), dummy_Ye, press(p.I), velx(p.I), vely(p.I),
velz(p.I), wlor, Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), Ex, Ey, Ez);

zvec_x(p.I) = wlor * pv.vel(0);
zvec_x(p.I) = wlor * pv.vel(0);
zvec_y(p.I) = wlor * pv.vel(1);
zvec_z(p.I) = wlor * pv.vel(2);

svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0);
svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1);
svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2);
svec_x(p.I) =
(pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(0);
svec_y(p.I) =
(pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(1);
svec_z(p.I) =
(pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(2);

// Write back cv
cv.scatter(dens(p.I), momx(p.I), momy(p.I), momz(p.I), tau(p.I), dummy_Ye,
Expand Down
134 changes: 126 additions & 8 deletions AsterX/src/con2prim_rpa.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,10 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
sm_metric3 g(sm_symt3l(glo(0, 0), glo(0, 1), glo(1, 1), glo(0, 2),
glo(1, 2), glo(2, 2)));

/* Get covariant metric used for BH interior fixes */
const smat<CCTK_REAL, 3> glow(
[&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); });

/* Calculate inverse of 3-metric */
const CCTK_REAL spatial_detg = calc_det(glo);
const CCTK_REAL sqrt_detg = sqrt(spatial_detg);
Expand All @@ -105,6 +109,66 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
{momx(p.I), momy(p.I), momz(p.I)},
{dBx(p.I), dBy(p.I), dBz(p.I)}};


// Modifying primitive seeds within BH interiors before C2Ps are called
// NOTE: By default, alp_thresh=0 so the if condition below is never
// triggered. One must be very careful when using this functionality and
// must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the parfile

if (alp(p.I) < alp_thresh) {
if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) {
pv_seeds.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of
// initial NS or disk
pv_seeds.eps = eps_BH;
pv_seeds.ye = Ye_atmo;
pv_seeds.press =
eos.at_rho_eps_ye(rho_BH, eps_BH, Ye_atmo).press();
if ( ((pv_seeds.vel(0) * pv_seeds.w_lor) > vwlim_BH) ||
((pv_seeds.vel(1)*pv_seeds.w_lor) > vwlim_BH) ||
((pv_seeds.vel(2)*pv_seeds.w_lor) > vwlim_BH) ) {
CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH);
CCTK_REAL vlim_BH = vwlim_BH / wlim_BH;
pv_seeds.vel(0) *= vlim_BH / pv_seeds.vel(0);
pv_seeds.vel(1) *= vlim_BH / pv_seeds.vel(1);
pv_seeds.vel(2) *= vlim_BH / pv_seeds.vel(2);
pv_seeds.w_lor = wlim_BH;
}
// cv.from_prim(pv_seeds, g);
// We do not save electric field, required by cv.from_prim(pv_seeds, g) in RPA.
// Thus, we recompute the CVs explicitly below

const vec<CCTK_REAL, 3> &v_up = pv_seeds.vel;
const vec<CCTK_REAL, 3> v_low = calc_contraction(glow, v_up);
/* Computing B_j */
const vec<CCTK_REAL, 3> &B_up = pv_seeds.B;
const vec<CCTK_REAL, 3> B_low = calc_contraction(glow, B_up);
/* Computing b^t : this is b^0 * alp */
const CCTK_REAL bst = pv_seeds.w_lor * calc_contraction(B_up, v_low);
/* Computing b_j */
const vec<CCTK_REAL, 3> b_low = B_low / pv_seeds.w_lor + bst * v_low;
/* Computing b^mu b_mu */
const CCTK_REAL bs2 =
(calc_contraction(B_up, B_low) + bst * bst) / (pv_seeds.w_lor * pv_seeds.w_lor);
// computing conserved from primitives
cv.dens = sqrt_detg * pv_seeds.rho * pv_seeds.w_lor;
cv.scon(0) = sqrt_detg * (pv_seeds.w_lor * pv_seeds.w_lor *
(pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * v_low(0) -
bst * b_low(0));
cv.scon(1) = sqrt_detg * (pv_seeds.w_lor * pv_seeds.w_lor *
(pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * v_low(1) -
bst * b_low(1));
cv.scon(2) = sqrt_detg * (pv_seeds.w_lor * pv_seeds.w_lor *
(pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * v_low(2) -
bst * b_low(2));
cv.tau = sqrt_detg * (pv_seeds.w_lor * pv_seeds.w_lor *
(pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) -
(pv_seeds.press + 0.5 * bs2) - bst * bst) - cv.dens;
cv.bcons = sqrt_detg * pv_seeds.B;
cv.tracer_ye = cv.dens * pv_seeds.ye;
}
}

// Calling C2P
cv2pv(pv, cv, g, rep);

// Handle incorrectable errors
Expand All @@ -115,8 +179,7 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
// need to fix pv to computed values like pv.rho instead of rho(p.I)
printf(
"WARNING: "
"C2P failed. Printing cons and saved prims before set to "
"atmo: \n"
"C2P failed. Printing cons and saved prims before set to atmo or BH interior fix: \n"
"cctk_iteration = %i \n "
"x, y, z = %26.16e, %26.16e, %26.16e \n "
"gxx, gxy, gxz, gyy, gyz, gzz = %f, %f, %f, %f, %f, %f \n "
Expand All @@ -138,12 +201,67 @@ 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
cv.bcons(0) = dBx(p.I);
cv.bcons(1) = dBy(p.I);
cv.bcons(2) = dBz(p.I);
pv.B = cv.bcons / sqrt_detg;
atmo.set(pv, cv, g);

if (alp(p.I) < alp_thresh) {
if ((pv.rho > rho_BH) || (pv.eps > eps_BH)) {
pv.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of
// initial NS or disk
pv.eps = eps_BH;
pv.ye = Ye_atmo;
pv.press =
eos.at_rho_eps_ye(rho_BH, eps_BH, Ye_atmo).press();
if (((pv.vel(0) * pv.w_lor) > vwlim_BH) ||
((pv.vel(1)*pv.w_lor) > vwlim_BH) ||
((pv.vel(2)*pv.w_lor) > vwlim_BH) ) {
CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH);
CCTK_REAL vlim_BH = vwlim_BH / wlim_BH;
pv.vel(0) *= vlim_BH / pv.vel(0);
pv.vel(1) *= vlim_BH / pv.vel(1);
pv.vel(2) *= vlim_BH / pv.vel(2);
pv.w_lor = wlim_BH;
}
// cv.from_prim(pv, g);
// We do not save electric field, required by cv.from_prim(pv, g) in RPA.
// Thus, we recompute the CVs explicitly below

const vec<CCTK_REAL, 3> &v_up = pv.vel;
const vec<CCTK_REAL, 3> v_low = calc_contraction(glow, v_up);
/* Computing B_j */
const vec<CCTK_REAL, 3> &B_up = pv.B;
const vec<CCTK_REAL, 3> B_low = calc_contraction(glow, B_up);
/* Computing b^t : this is b^0 * alp */
const CCTK_REAL bst = pv.w_lor * calc_contraction(B_up, v_low);
/* Computing b_j */
const vec<CCTK_REAL, 3> b_low = B_low / pv.w_lor + bst * v_low;
/* Computing b^mu b_mu */
const CCTK_REAL bs2 =
(calc_contraction(B_up, B_low) + bst * bst) / (pv.w_lor * pv.w_lor);
// computing conserved from primitives
cv.dens = sqrt_detg * pv.rho * pv.w_lor;
cv.scon(0) = sqrt_detg * (pv.w_lor * pv.w_lor *
(pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(0) -
bst * b_low(0));
cv.scon(1) = sqrt_detg * (pv.w_lor * pv.w_lor *
(pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(1) -
bst * b_low(1));
cv.scon(2) = sqrt_detg * (pv.w_lor * pv.w_lor *
(pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(2) -
bst * b_low(2));
cv.tau = sqrt_detg * (pv.w_lor * pv.w_lor *
(pv.rho * (1.0 + pv.eps) + pv.press + bs2) -
(pv.press + 0.5 * bs2) - bst * bst) - cv.dens;
cv.bcons = sqrt_detg * pv.B;
cv.tracer_ye = cv.dens * pv.ye;
}
}
else {
// set to atmo
cv.bcons(0) = dBx(p.I);
cv.bcons(1) = dBy(p.I);
cv.bcons(2) = dBz(p.I);
pv.B = cv.bcons / sqrt_detg;
atmo.set(pv, cv, g);
}
}

/* set flag to success */
Expand Down

0 comments on commit 9e28ade

Please sign in to comment.