diff --git a/AsterSeeds/interface.ccl b/AsterSeeds/interface.ccl index 925cc465..5c8d0b82 100644 --- a/AsterSeeds/interface.ccl +++ b/AsterSeeds/interface.ccl @@ -1,7 +1,7 @@ # Interface definition for thorn AsterSeeds IMPLEMENTS: AsterSeeds -INHERITS: HydroBaseX AsterX +INHERITS: HydroBaseX AsterX VolumeIntegrals_GRMHDX USES INCLUDE HEADER: loop_device.hxx USES INCLUDE HEADER: setup_eos.hxx @@ -14,3 +14,31 @@ CCTK_REAL Avec_cent TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' Avec_z_cent } "Cell-centered vector potential components" +CCTK_REAL vel_NS1 TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=3 +{ + vel_NS1 +} "Helper for dynamic seeding of external dipole field. Records velocity of NS CoM" + +CCTK_REAL vel_NS2 TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=3 +{ + vel_NS2 +} "Helper for dynamic seeding of external dipole field. Records velocity of NS CoM" + +CCTK_INT SeedNow type=SCALAR tags='checkpoint="yes"' +CCTK_INT DoneSeeding type=SCALAR tags='checkpoint="yes"' + +CCTK_INT FUNCTION DriverInterpolate( + CCTK_POINTER_TO_CONST IN cctkGH, + CCTK_INT IN N_dims, + CCTK_INT IN local_interp_handle, + CCTK_INT IN param_table_handle, + CCTK_INT IN coord_system_handle, + CCTK_INT IN N_interp_points, + CCTK_INT IN interp_coords_type, + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, + CCTK_INT IN N_input_arrays, + CCTK_INT ARRAY IN input_array_indices, + CCTK_INT IN N_output_arrays, + CCTK_INT ARRAY IN output_array_types, + CCTK_POINTER ARRAY IN output_arrays) +REQUIRES FUNCTION DriverInterpolate diff --git a/AsterSeeds/param.ccl b/AsterSeeds/param.ccl index 2b3cc001..5dfdf4df 100644 --- a/AsterSeeds/param.ccl +++ b/AsterSeeds/param.ccl @@ -165,6 +165,7 @@ KEYWORD Afield_config "Definition of the initial vector potential" "none" :: "Nothing set here" "internal dipole" :: "Dipole field according to Ciolfi+2017" "external dipole" :: "Dipole field according to Moesta+2020" + "external dipole UIUC" :: "Dipole field according to Paschalidis+2013" } "none" # parameters for internal dipolar magnetic field @@ -218,6 +219,65 @@ REAL dipole_z[2] "z-coordinate of the dipole center" STEERABLE=ALWAYS *:* :: "Anything" } 0.0 +# Initial low density magnetosphere parameters + +BOOLEAN initial_beta "Pad pressure to set initial plasma beta floor?" STEERABLE=ALWAYS +{ +} "no" + +REAL initial_beta_min "Initial plasma beta floor" STEERABLE=ALWAYS +{ + 0:* :: "Positive" +} 0.0 + +REAL initial_beta_rhocut "Pad pressure below this density" STEERABLE=ALWAYS +{ + 0:* :: "Positive" +} 0.0 + +BOOLEAN set_coorbiting_vel "After padding pressure, set coorbiting velocity? This feature is not well tested, use with caution." STEERABLE=ALWAYS +{ +} "no" + +REAL radius_NS1 "Approximate radius of initial left NS, used only for tapering coorbiting vel." STEERABLE=ALWAYS +{ + 0:* :: "Positive" +} 9.0 + +REAL radius_NS2 "Approximate radius of initial right NS, used only for tapering coorbiting vel." STEERABLE=ALWAYS +{ + 0:* :: "Positive" +} 9.0 + +REAL vtol "Overwrite velocity with coorbiting velocity if all |v^i| are below tolerance" STEERABLE=ALWAYS +{ + 0:* :: "Positive" +} 0.0 + +# BNS inspiral initialization parameters + +BOOLEAN use_separation "Seeding time determined by separation of stars from VolumeIntegrals_GRMHDX" STEERABLE=ALWAYS +{ +} "yes" + +BOOLEAN use_time "Seeding time determined by evolution time" STEERABLE=ALWAYS +{ +} "no" + +REAL seeding_separation "Coordinate separation of stars below which magnetic field is seeded" STEERABLE=ALWAYS +{ + (0.0:*) :: "" +} 5.0 + +REAL seeding_time "Time [Cactus units] after which magnetic field is seeded" STEERABLE=ALWAYS +{ + [0.0:*) :: "" +} 0.0 + +INT seed_every "Check if conditions for seeding are met every iterations" STEERABLE=ALWAYS +{ + 0:* :: "0 means dont seed in CCTK_POSTSTEP" +} 0 # --------------------------------------- @@ -311,3 +371,7 @@ USES KEYWORD evolution_eos SHARES: Con2PrimFactory USES CCTK_REAL Ye_atmo +USES CCTK_REAL rho_abs_min +USES CCTK_REAL atmo_tol +USES CCTK_REAL r_atmo +USES CCTK_REAL n_rho_atmo diff --git a/AsterSeeds/schedule.ccl b/AsterSeeds/schedule.ccl index 0d02cc5f..6638d48d 100644 --- a/AsterSeeds/schedule.ccl +++ b/AsterSeeds/schedule.ccl @@ -87,22 +87,148 @@ if (CCTK_Equals(test_type, "3DTest")) { if (CCTK_Equals(test_case, "magBNS")) { STORAGE: Avec_cent - SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS AT initial IN ODESolvers_Initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup { - LANG: C - READS: HydroBaseX::press(everywhere) - WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + SCHEDULE AsterSeeds_Initialize_Seeding_Flags AT CCTK_INITIAL + { + LANG: C + OPTIONS: GLOBAL + WRITES: SeedNow(everywhere), DoneSeeding(everywhere) + WRITES: vel_NS1(everywhere) vel_NS2(everywhere) + } "Initialize seeding variables." + + if (seed_every != 0) { + SCHEDULE AsterSeeds_InitializeAvectoZero AT initial IN ODESolvers_Initial BEFORE AsterX_InitialGroup + { + LANG: C + WRITES: AsterX::Avec_x(everywhere) AsterX::Avec_y(everywhere) AsterX::Avec_z(everywhere) + } "Set up zero vector potential" + + SCHEDULE GROUP Seed_During_Evolution at CCTK_POSTSTEP + { + } "Routines for seeding magnetic field during evolution." + + SCHEDULE AsterSeeds_Set_Seeding_Flags in Seed_During_Evolution + { + LANG: C + OPTIONS: GLOBAL + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + READS: DoneSeeding(everywhere) + WRITES: SeedNow(everywhere) + } "Evaluate whether to seed now or not." + + SCHEDULE GROUP Seeding_Tasks IN Seed_During_Evolution IF AsterSeeds::SeedNow AFTER AsterSeeds_Set_Seeding_Flags + { + } "If conditions are met, seed magnetic field." + + SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS IN Seeding_Tasks + { + LANG: C + READS: HydroBaseX::press(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + WRITES: SeedNow(everywhere), DoneSeeding(everywhere) + } "Set up initial conditions for the cell-centered vector potential" + + SCHEDULE AsterSeeds_InitializeStagAvec_BNS IN Seeding_Tasks AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup + { + LANG: C + READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior) + SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z + } "Set up initial conditions for the vector potential" + + SCHEDULE GROUP AsterX_Con2PrimGroup IN Seeding_Tasks AFTER AsterSeeds_InitializeStagAvec_BNS + { + } "Calc B from A, run Con2Prim." + + if (initial_beta) { + if (set_coorbiting_vel) { + SCHEDULE AsterSeeds_InterpolateNSVelocity IN Seeding_Tasks BEFORE AsterSeeds_SetInitialBetaFloor + { + LANG: C + OPTIONS: GLOBAL + READS: HydroBaseX::vel(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + WRITES: vel_NS1(everywhere) vel_NS2(everywhere) + } "Get velocity of NS CoMs" + } + + SCHEDULE AsterSeeds_SetInitialBetaFloor IN Seeding_Tasks AFTER AsterX_Con2PrimGroup + { + LANG: C + READS: HydroBaseX::vel(everywhere) HydroBaseX::Bvec(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + READS: ADMBaseX::metric(everywhere) + READS: vel_NS1(everywhere) vel_NS2(everywhere) + WRITES: HydrobaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::entropy(everywhere) HydrobaseX::vel(everywhere) + } "Enforce initial minimum plasma beta" + + SCHEDULE AsterX_Prim2Con_Initial IN Seeding_Tasks AFTER AsterSeeds_SetInitialBetaFloor + { + LANG: C + READS: ADMBaseX::metric(interior) + READS: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior) + READS: HydroBaseX::entropy(interior) HydroBaseX::Ye(interior) + WRITES: AsterX::cons_vector(interior) AsterX::dB(interior) + SYNC: AsterX::cons_vector AsterX::dB + } "Update cons after setting initial beta" + } + + SCHEDULE AsterX_Tmunu IN Seeding_Tasks AFTER AsterX_Con2PrimGroup AFTER AsterX_ApplyOuterBCOnPrim + { + LANG: C + READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere) + READS: HydroBaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) + READS: HydroBaseX::vel(everywhere) + READS: AsterX::zvec_x(everywhere), AsterX::zvec_y(everywhere), AsterX::zvec_z(everywhere) + READS: AsterX::svec_x(everywhere), AsterX::svec_y(everywhere), AsterX::svec_z(everywhere) + READS: HydroBaseX::Bvec(everywhere) + READS: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior) + WRITES: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior) + SYNC: TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij + } "Compute the energy-momentum tensor" + } - "Set up initial conditions for the cell-centered vector potential" + else { + SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS AT initial IN ODESolvers_Initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup { + LANG: C + READS: HydroBaseX::press(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + WRITES: SeedNow(everywhere), DoneSeeding(everywhere) + } "Set up initial conditions for the cell-centered vector potential" - SCHEDULE AsterSeeds_InitializeStagAvec_BNS AT initial IN ODESolvers_Initial AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup { - LANG: C - READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) - WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior) - SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z + SCHEDULE AsterSeeds_InitializeStagAvec_BNS AT initial IN ODESolvers_Initial AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup + { + LANG: C + READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior) + SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z + } "Set up initial conditions for the vector potential" + + if (initial_beta) { + if (set_coorbiting_vel) { + SCHEDULE AsterSeeds_InterpolateNSVelocity IN Seeding_Tasks BEFORE AsterSeeds_SetInitialBetaFloor + { + LANG: C + OPTIONS: GLOBAL + READS: HydroBaseX::vel(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + WRITES: vel_NS1(everywhere) vel_NS2(everywhere) + } "Get velocity of NS CoMs" + } + + SCHEDULE AsterSeeds_SetInitialBetaFloor AT initial IN AsterX_InitialGroup AFTER AsterX_ComputeBFromdB BEFORE AsterX_CheckPrims + { + LANG: C + READS: HydroBaseX::vel(everywhere) HydroBaseX::Bvec(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere) + READS: ADMBaseX::metric(everywhere) + READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere) + READS: vel_NS1(everywhere) vel_NS2(everywhere) + WRITES: HydrobaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::entropy(everywhere) HydrobaseX::vel(everywhere) + } "Enforce initial minimum plasma beta" + } } - "Set up initial conditions for the vector potential" } - } #Initial conditions for TabEOS tests diff --git a/AsterSeeds/src/initial_betafloor.cxx b/AsterSeeds/src/initial_betafloor.cxx new file mode 100644 index 00000000..4cc8e368 --- /dev/null +++ b/AsterSeeds/src/initial_betafloor.cxx @@ -0,0 +1,189 @@ +#include + +#include +#include +#include + +#include + +#include "util_Table.h" +#include "seeds_utils.hxx" +#include "setup_eos.hxx" + +namespace AsterSeeds { +using namespace std; +using namespace Loop; +using namespace AsterUtils; +using namespace EOSX; + +extern "C" void AsterSeeds_InterpolateNSVelocity(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InterpolateNSVelocity; + DECLARE_CCTK_PARAMETERS; + + // Get NS velocities + const int nPoints = 2; + const int nInputArrays = 3; + CCTK_REAL nsx[nPoints] = {*comx1, *comx2}; + CCTK_REAL nsy[nPoints] = {*comy1, *comy2}; + CCTK_REAL nsz[nPoints] = {*comz1, *comz2}; + const void *interp_coords[nInputArrays] = {(const void *)nsx, (const void *)nsy, + (const void *)nsz}; + const CCTK_INT inputArrayIndices[nInputArrays] = { + CCTK_VarIndex("HydroBaseX::velx"), CCTK_VarIndex("HydroBaseX::vely"), CCTK_VarIndex("HydroBaseX::velz")}; + CCTK_REAL nsvx[nPoints], nsvy[nPoints], nsvz[nPoints]; + CCTK_POINTER outputArrays[nInputArrays] = {(void*)nsvx, (void*)nsvy, (void*)nsvz}; + + // DriverInterpolate arguments that aren't currently used + const int coordSystemHandle = 0; + const CCTK_INT interpCoordsTypeCode = 0; + const CCTK_INT outputArrayTypes[nInputArrays] = {0, 0, 0}; + + const int interpHandle = CCTK_InterpHandle("CarpetX"); + if (interpHandle < 0) { + CCTK_WARN(CCTK_WARN_ALERT, "Can't get interpolation handle"); + return; + } + + // Create parameter table for interpolation + const int paramTableHandle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); + if (paramTableHandle < 0) { + CCTK_VERROR("Can't create parameter table: %d", paramTableHandle); + } + + // Set interpolation order in the parameter table + int ierr = Util_TableSetInt(paramTableHandle, 1, "order"); + if (ierr < 0) { + CCTK_VERROR("Can't set order in parameter table: %d", ierr); + } + + // Perform the interpolation + ierr = DriverInterpolate(cctkGH, 3, interpHandle, paramTableHandle, + coordSystemHandle, nPoints, interpCoordsTypeCode, + interp_coords, nInputArrays, inputArrayIndices, + nInputArrays, outputArrayTypes, outputArrays); + + CCTK_VINFO("Interpolated (%g, %g, %g) as NS1 velocity", nsvx[0], nsvy[0], nsvz[0]); + CCTK_VINFO("Interpolated (%g, %g, %g) as NS2 velocity", nsvx[1], nsvy[1], nsvz[1]); + + vel_NS1[0] = nsvx[0]; + vel_NS1[1] = nsvy[0]; + vel_NS1[2] = nsvz[0]; + vel_NS2[0] = nsvx[1]; + vel_NS2[1] = nsvy[1]; + vel_NS2[2] = nsvz[1]; + + return; +} + +extern "C" void AsterSeeds_SetInitialBetaFloor(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_SetInitialBetaFloor; + DECLARE_CCTK_PARAMETERS; + + auto eos_3p_tab3d = global_eos_3p_tab3d; + if (not CCTK_EQUALS(evolution_eos, "Tabulated3d")) { + CCTK_VERROR("Invalid evolution EOS type '%s'. Please, set " + "EOSX::evolution_eos = \"Tabulated3d\" in your parameter file.", + evolution_eos); + } + + const smat, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz}; + + CCTK_VINFO("Using (%g, %g, %g) as NS1 velocity", vel_NS1[0], vel_NS1[1], vel_NS1[2]); + CCTK_VINFO("Using (%g, %g, %g) as NS2 velocity", vel_NS2[0], vel_NS2[1], vel_NS2[2]); + + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + + const CCTK_REAL pressL = press(p.I); + const CCTK_REAL rhoL = rho(p.I); + const CCTK_REAL tempL = temperature(p.I); + const CCTK_REAL YeL = Ye(p.I); + const CCTK_REAL epsL = eps(p.I); + const CCTK_REAL entL = entropy(p.I); + + // Compute b^2 + /* Get covariant metric */ + const smat glo( + [&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); }); + + vec B_up{Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)}; + vec B_low = calc_contraction(glo, B_up); + + vec v_up{velx(p.I), vely(p.I), velz(p.I)}; + vec v_low = calc_contraction(glo, v_up); + + const CCTK_REAL wlor = calc_wlorentz(v_up, v_low); + const CCTK_REAL alp_b0 = wlor * calc_contraction(B_up, v_low); + + const CCTK_REAL B2 = calc_contraction(B_up, B_low); + const CCTK_REAL bsq = ( B2 + alp_b0 * alp_b0 ) / ( wlor*wlor ); + + // Increase P if necessary + CCTK_REAL press_lim = initial_beta_min * bsq * 0.5; + + if ((pressL >= press_lim) || (rhoL > initial_beta_rhocut)) { + press(p.I) = pressL; + rho(p.I) = rhoL; + eps(p.I) = epsL; + entropy(p.I) = entL; + } + else { + // Recalculate primitives + press(p.I) = press_lim; + rho(p.I) = eos_3p_tab3d->rho_from_press_temp_ye(press_lim, tempL, YeL); + eps(p.I) = eos_3p_tab3d->eps_from_rho_temp_ye(rho(p.I), tempL, YeL); + entropy(p.I) = eos_3p_tab3d->entropy_from_rho_temp_ye(rho(p.I), tempL, YeL); + } + + // TODO: The coorbiting velocity feature is not well tested. Use with caution. + if (set_coorbiting_vel) { + const CCTK_REAL vxL = velx(p.I); + const CCTK_REAL vyL = vely(p.I); + const CCTK_REAL vzL = velz(p.I); + if ((abs(vxL) < vtol) && (abs(vyL) < vtol) && (abs(vzL) < vtol)) { + // For star 1 at minus side + const CCTK_REAL x_local_s1 = p.x - *comx1; + const CCTK_REAL y_local_s1 = p.y - *comy1; + const CCTK_REAL cylrad2_s1 = + x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1; + // For star 2 at minus side + const CCTK_REAL x_local_s2 = p.x - *comx2; + const CCTK_REAL y_local_s2 = p.y - *comy2; + const CCTK_REAL cylrad2_s2 = + x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2; + + if (cylrad2_s1 < cylrad2_s2) { + if (cylrad2_s1 < pow(3.0 * radius_NS1, 2.0)) { + velx(p.I) = vel_NS1[0]; + vely(p.I) = vel_NS1[1]; + velz(p.I) = vel_NS1[2]; + } else { + CCTK_REAL rfac = pow(3.0 * radius_NS1, 4.0) / pow(cylrad2_s1, 2.0); + velx(p.I) = rfac * vel_NS1[0]; + vely(p.I) = rfac * vel_NS1[1]; + velz(p.I) = rfac * vel_NS1[2]; + } + } else { + if (cylrad2_s2 < pow(3.0 * radius_NS2, 2.0)) { + velx(p.I) = vel_NS2[0]; + vely(p.I) = vel_NS2[1]; + velz(p.I) = vel_NS2[2]; + } else { + CCTK_REAL rfac = pow(3.0 * radius_NS2, 4.0) / pow(cylrad2_s2, 2.0); + velx(p.I) = rfac * vel_NS2[0]; + vely(p.I) = rfac * vel_NS2[1]; + velz(p.I) = rfac * vel_NS2[2]; + } + } + } else { + velx(p.I) = vxL; + vely(p.I) = vyL; + velz(p.I) = vzL; + } + } // if set coorbiting velocity + + }); +} + +} // namespace AsterSeeds diff --git a/AsterSeeds/src/magBNS.cxx b/AsterSeeds/src/magBNS.cxx index 899fbae7..d226456b 100644 --- a/AsterSeeds/src/magBNS.cxx +++ b/AsterSeeds/src/magBNS.cxx @@ -13,10 +13,67 @@ using namespace std; using namespace Loop; using namespace AsterUtils; +extern "C" void AsterSeeds_Initialize_Seeding_Flags(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_Initialize_Seeding_Flags; + DECLARE_CCTK_PARAMETERS; + + *SeedNow = 0; + *DoneSeeding = 0; + for (int ii=0; ii<3; ii++) { + vel_NS1[ii] = 0.0; + vel_NS2[ii] = 0.0; + } + return; +} + +extern "C" void AsterSeeds_Set_Seeding_Flags(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_Set_Seeding_Flags; + DECLARE_CCTK_PARAMETERS; + + bool seed_trigger = false; + if (use_time && (cctk_time > seeding_time)) + seed_trigger = true; + if (use_separation) { + const CCTK_REAL ns_Dx = abs(*comx1 - *comx2); + const CCTK_REAL ns_Dy = abs(*comy1 - *comy2); + const CCTK_REAL ns_sep = sqrt(ns_Dx * ns_Dx + ns_Dy * ns_Dy); + if (ns_sep > seeding_separation) + seed_trigger = true; + } + + if ((!*DoneSeeding) && (cctk_iteration % seed_every == 0) && seed_trigger) { + CCTK_VINFO("Seeding Initial Avec..."); + *SeedNow = 1; + } else { + *SeedNow = 0; + } + + return; +} + extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeCenteredAvec_BNS; DECLARE_CCTK_PARAMETERS; + // Set origin according to parameter or NS CoMs + CCTK_REAL x01, y01, z01, x02, y02, z02; + if (seed_every != 0) { + x01 = *comx1; + y01 = *comy1; + z01 = *comz1; + x02 = *comx2; + y02 = *comy2; + z02 = *comz2; + } + else { + x01 = dipole_x[0]; + y01 = dipole_y[0]; + z01 = dipole_z[0]; + x02 = dipole_x[1]; + y02 = dipole_y[1]; + z02 = dipole_z[1]; + } + if (CCTK_EQUALS(Afield_config, "internal dipole")) { /* computing cell centered vector potential components */ @@ -29,13 +86,13 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { // For star 1 at minus side if (p.x < 0) { - x_local = p.x - dipole_x[0]; - y_local = p.y - dipole_y[0]; + x_local = p.x - x01; + y_local = p.y - y01; } - // For star 2 at minus side + // For star 2 at plus side else { - x_local = p.x - dipole_x[1]; - y_local = p.y - dipole_y[1]; + x_local = p.x - x02; + y_local = p.y - y02; } CCTK_REAL Pcut = press_max * press_cut; @@ -53,9 +110,9 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { grid.nghostzones, [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { // For star 1 at minus side - CCTK_REAL x_local_s1 = p.x - dipole_x[0]; - CCTK_REAL y_local_s1 = p.y - dipole_y[0]; - CCTK_REAL z_local_s1 = p.z - dipole_z[0]; + CCTK_REAL x_local_s1 = p.x - x01; + CCTK_REAL y_local_s1 = p.y - y01; + CCTK_REAL z_local_s1 = p.z - z01; CCTK_REAL cylrad2_s1 = x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1; CCTK_REAL rsph_s1 = @@ -67,9 +124,9 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { sqrt(cylrad2_s1 + 1.0e-16); // For star 2 at minus side - CCTK_REAL x_local_s2 = p.x - dipole_x[1]; - CCTK_REAL y_local_s2 = p.y - dipole_y[1]; - CCTK_REAL z_local_s2 = p.z - dipole_z[1]; + CCTK_REAL x_local_s2 = p.x - x02; + CCTK_REAL y_local_s2 = p.y - y02; + CCTK_REAL z_local_s2 = p.z - z02; CCTK_REAL cylrad2_s2 = x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2; CCTK_REAL rsph_s2 = @@ -87,9 +144,57 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { Avec_z_cent(p.I) = 0.0; }); + } else if (CCTK_EQUALS(Afield_config, "external dipole UIUC")) { + + /* computing cell centered vector potential components */ + grid.loop_all<1, 1, 1>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + + const CCTK_REAL pi = 2 * acos(0.0); + + // For star 1 at minus side + CCTK_REAL x_local_s1 = p.x - x01; + CCTK_REAL y_local_s1 = p.y - y01; + CCTK_REAL z_local_s1 = p.z - z01; + CCTK_REAL cylrad2_s1 = + x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1; + CCTK_REAL sphrad2_s1 = + x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1 + + z_local_s1 * z_local_s1; + CCTK_REAL r02 = r0 * r0; + + // See e.g. Ruiz+ 2020, Eq. 1. Here, B0 is I0 from Eq. 1, and we have preemptively canceled out the factor + // of cylrad2. + CCTK_REAL Aphi_local_s1 = + pi * B0 * r02 / pow(r02 + sphrad2_s1, 1.5) + * (1.0 + (15.0 * r02 * (r02 + cylrad2_s1) / (8.0 * pow(r02 + sphrad2_s1, 2.0)))); + + // For star 2 at minus side + CCTK_REAL x_local_s2 = p.x - x02; + CCTK_REAL y_local_s2 = p.y - y02; + CCTK_REAL z_local_s2 = p.z - z02; + CCTK_REAL cylrad2_s2 = + x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2; + CCTK_REAL sphrad2_s2 = + x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2 + + z_local_s2 * z_local_s2; + + CCTK_REAL Aphi_local_s2 = + pi * B0 * r02 / pow(r02 + sphrad2_s2, 1.5) + * (1.0 + (15.0 * r02 * (r02 + cylrad2_s2) / (8.0 * pow(r02 + sphrad2_s2, 2.0)))); + + Avec_x_cent(p.I) = + -(y_local_s1 * Aphi_local_s1 + y_local_s2 * Aphi_local_s2); + Avec_y_cent(p.I) = + x_local_s1 * Aphi_local_s1 + x_local_s2 * Aphi_local_s2; + Avec_z_cent(p.I) = 0.0; + }); } else { CCTK_ERROR("Vector potential configuration not defined"); } + + *DoneSeeding = 1; } extern "C" void AsterSeeds_InitializeStagAvec_BNS(CCTK_ARGUMENTS) { @@ -115,4 +220,27 @@ extern "C" void AsterSeeds_InitializeStagAvec_BNS(CCTK_ARGUMENTS) { }); } +extern "C" void AsterSeeds_InitializeAvectoZero(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeAvectoZero; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all<1, 0, 0>(grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_x(p.I) = 0.0; + }); + + grid.loop_all<0, 1, 0>(grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_y(p.I) = 0.0; + }); + + grid.loop_all<0, 0, 1>(grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_z(p.I) = 0.0; + }); +} + } // namespace AsterSeeds diff --git a/AsterSeeds/src/make.code.defn b/AsterSeeds/src/make.code.defn index a9bd4c0d..96819da3 100644 --- a/AsterSeeds/src/make.code.defn +++ b/AsterSeeds/src/make.code.defn @@ -8,6 +8,7 @@ SRCS = \ TabEOS_tests.cxx \ atmosphere.cxx \ bbhcloud.cxx \ + initial_betafloor.cxx \ magTOV.cxx \ magBNS.cxx \ SetTempEntropyYe.cxx diff --git a/AsterX/interface.ccl b/AsterX/interface.ccl index 10b95dcb..634cfe1a 100644 --- a/AsterX/interface.ccl +++ b/AsterX/interface.ccl @@ -17,12 +17,12 @@ PUBLIC: # Note that the prolongation method and order are hard-coded and will # supersede parfile choices -CCTK_REAL cons_vector TYPE=gf CENTERING={ccc} TAGS='rhs="cons_vectorrhs" dependents="HydroBaseX::rho HydroBaseX::eps HydroBaseX::press HydroBaseX::temperature HydroBaseX::Ye HydroBaseX::Bvec HydroBaseX::vel HydroBaseX::entropy TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij AsterX::flux_x AsterX::flux_y AsterX::flux_z" prolongation_type="minmod" prolongation_order=1' "Conserved quantities during AsterX evolution" +CCTK_REAL cons_vector TYPE=gf CENTERING={ccc} TAGS='rhs="cons_vector_rhs" dependents="HydroBaseX::rho HydroBaseX::eps HydroBaseX::press HydroBaseX::temperature HydroBaseX::Ye HydroBaseX::Bvec HydroBaseX::vel HydroBaseX::entropy TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij AsterX::flux_x AsterX::flux_y AsterX::flux_z" prolongation_type="minmod" prolongation_order=1' "Conserved quantities during AsterX evolution" { dens, momx, momy, momz, tau, DYe, DEnt } "Conserved density, momenta, internal energy, electron fraction and entropy" -CCTK_REAL cons_vectorrhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "RHSs of conserved variables" +CCTK_REAL cons_vector_rhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "RHSs of conserved variables" { densrhs, momxrhs, momyrhs, momzrhs, taurhs, DYe_rhs, DEntrhs } "Conserved density, momenta, internal energy, electron fraction and entropy RHS" @@ -140,6 +140,7 @@ CCTK_REAL diagnostics TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' w_lorentz B_norm b2small + volform } "Diagnostic grid functions" # saved_prims are computed only once at the beginning of the simulation @@ -199,7 +200,7 @@ CCTK_REAL a_zface TYPE=gf CENTERING={ccv} TAGS='checkpoint="no"' # For Subcycling -CCTK_REAL cons_vectorold TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' +CCTK_REAL cons_vector_old TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' { densold, momxold, momyold, momzold, tauold, DYeold, DEntold } "Conserved vector RHS" @@ -208,7 +209,7 @@ CCTK_REAL Avec_y_old TYPE=gf CENTERING={vcv} TAGS='checkpoint="no" prolongation_ CCTK_REAL Avec_z_old TYPE=gf CENTERING={vvc} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "z-component of vector potential RHS" CCTK_REAL Psi_old TYPE=gf CENTERING={vvv} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "t-component of vector potential RHS" -CCTK_REAL cons_vectork1 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' +CCTK_REAL cons_vector_k1 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' { densk1, momxk1, momyk1, momzk1, tauk1, DYek1, DEntk1 } "Conserved vector RHS" @@ -217,7 +218,7 @@ CCTK_REAL Avec_y_k1 TYPE=gf CENTERING={vcv} TAGS='checkpoint="no" prolongation_t CCTK_REAL Avec_z_k1 TYPE=gf CENTERING={vvc} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "z-component of vector potential RHS" CCTK_REAL Psi_k1 TYPE=gf CENTERING={vvv} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "t-component of vector potential RHS" -CCTK_REAL cons_vectork2 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' +CCTK_REAL cons_vector_k2 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' { densk2, momxk2, momyk2, momzk2, tauk2, DYek2, DEntk2 } "Conserved vector RHS" @@ -226,7 +227,7 @@ CCTK_REAL Avec_y_k2 TYPE=gf CENTERING={vcv} TAGS='checkpoint="no" prolongation_t CCTK_REAL Avec_z_k2 TYPE=gf CENTERING={vvc} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "z-component of vector potential RHS" CCTK_REAL Psi_k2 TYPE=gf CENTERING={vvv} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "t-component of vector potential RHS" -CCTK_REAL cons_vectork3 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' +CCTK_REAL cons_vector_k3 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' { densk3, momxk3, momyk3, momzk3, tauk3, DYek3, DEntk3 } "Conserved vector RHS" @@ -235,7 +236,7 @@ CCTK_REAL Avec_y_k3 TYPE=gf CENTERING={vcv} TAGS='checkpoint="no" prolongation_t CCTK_REAL Avec_z_k3 TYPE=gf CENTERING={vvc} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "z-component of vector potential RHS" CCTK_REAL Psi_k3 TYPE=gf CENTERING={vvv} TAGS='checkpoint="no" prolongation_type="hermite" prolongation_order=5' "t-component of vector potential RHS" -CCTK_REAL cons_vectork4 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' +CCTK_REAL cons_vector_k4 TYPE=gf CENTERING={ccc} TAGS='checkpoint="no" prolongation_type="minmod" prolongation_order=1' { densk4, momxk4, momyk4, momzk4, tauk4, DYek4, DEntk4 } "Conserved vector RHS" diff --git a/AsterX/schedule.ccl b/AsterX/schedule.ccl index 5c053459..8bf61f89 100644 --- a/AsterX/schedule.ccl +++ b/AsterX/schedule.ccl @@ -2,7 +2,7 @@ STORAGE: cons_vector dB Psi HydroBaseX::Bvec dBx_stag dBy_stag dBz_stag STORAGE: flux_x flux_y flux_z -STORAGE: cons_vectorrhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs +STORAGE: cons_vector_rhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs STORAGE: ADMBaseX::metric ADMBaseX::lapse ADMBaseX::shift ADMBaseX::curv STORAGE: G Fbeta Fx_stag Fy_stag Fz_stag STORAGE: TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij @@ -10,11 +10,11 @@ STORAGE: Ex Ey Ez STORAGE: con2prim_flag diagnostics if(use_subcycling_wip) { - STORAGE: cons_vectorold Avec_x_old Avec_y_old Avec_z_old Psi_old - STORAGE: cons_vectork1 Avec_x_k1 Avec_y_k1 Avec_z_k1 Psi_k1 - STORAGE: cons_vectork2 Avec_x_k2 Avec_y_k2 Avec_z_k2 Psi_k2 - STORAGE: cons_vectork3 Avec_x_k3 Avec_y_k3 Avec_z_k3 Psi_k3 - STORAGE: cons_vectork4 Avec_x_k4 Avec_y_k4 Avec_z_k4 Psi_k4 + STORAGE: cons_vector_old Avec_x_old Avec_y_old Avec_z_old Psi_old + STORAGE: cons_vector_k1 Avec_x_k1 Avec_y_k1 Avec_z_k1 Psi_k1 + STORAGE: cons_vector_k2 Avec_x_k2 Avec_y_k2 Avec_z_k2 Psi_k2 + STORAGE: cons_vector_k3 Avec_x_k3 Avec_y_k3 Avec_z_k3 Psi_k3 + STORAGE: cons_vector_k4 Avec_x_k4 Avec_y_k4 Avec_z_k4 Psi_k4 SCHEDULE GROUP AsterX_Con2PrimGroup IN ODESolvers_PostSubStep BEFORE TmunuBaseX_SetTmunuVars AFTER ADMBaseX_SetADMVars { @@ -351,7 +351,7 @@ SCHEDULE AsterX_SourceTerms IN AsterX_RHSGroup AFTER AsterX_Fluxes READS: zvec_x(everywhere), zvec_y(everywhere), zvec_z(everywhere) READS: svec_x(everywhere), svec_y(everywhere), svec_z(everywhere) READS: HydroBaseX::Bvec(everywhere) - WRITES: cons_vectorrhs(interior) + WRITES: cons_vector_rhs(interior) } "Calculate the source terms and compute the RHS of the hydro equations" SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms @@ -359,7 +359,7 @@ SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms LANG: C READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere) READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere) - READS: cons_vectorrhs(interior) + READS: cons_vector_rhs(interior) READS: Psi(everywhere) READS: Ex(interior) Ey(interior) Ez(interior) READS: G(everywhere) @@ -367,7 +367,7 @@ SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms READS: Fx_stag(everywhere) Fy_stag(everywhere) Fz_stag(everywhere) READS: theta_x(everywhere) theta_y(everywhere) theta_z(everywhere) READS: LOflag(everywhere) - WRITES: cons_vectorrhs(interior) + WRITES: cons_vector_rhs(interior) WRITES: theta_tot(interior) WRITES: Avec_x_rhs(interior) Avec_y_rhs(interior) Avec_z_rhs(interior) Psi_rhs(interior) } "Update the RHS of the hydro equations with the flux contributions" diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 886e55be..b74f4968 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -501,6 +501,7 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType *eos_1p, w_lorentz(p.I) = wlor; B_norm(p.I) = sqrt(B2); b2small(p.I) = bsq; + volform(p.I) = sqrt_detg; }; cctk_grid.loop_all_device<1, 1, 1>(grid.nghostzones, c2p_impl); diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx index e92a4081..624e2290 100644 --- a/Con2PrimFactory/src/c2p.hxx +++ b/Con2PrimFactory/src/c2p.hxx @@ -254,17 +254,17 @@ c2p::prims_floors_and_ceilings(const EOSType *eos_3p, prim_vars &pv, rep.adjust_cons = true; if (use_temp) { - // Revert pressure change above for tabulated EOS - // TODO: Add functionality to support inv_beta limit w - // tabulated EOS + // Recompute T from adjusted rho, P + pv.temperature = eos_3p->temp_from_rho_press_ye(pv.rho, pv.press, pv.Ye); pv.eps = eos_3p->eps_from_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); - pv.press = eos_3p->press_from_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); - } else { + pv.entropy = eos_3p->entropy_from_rho_temp_ye(pv.rho, pv.temperature, pv.Ye); + } + else { pv.eps = eos_3p->eps_from_rho_press_ye(pv.rho, pv.press, pv.Ye); pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); } - pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); mag_ceiling = false; // Drift floors from https://arxiv.org/pdf/1611.09365 diff --git a/EOSX/src/eos_3p_hybrid/eos_3p_hybrid.hxx b/EOSX/src/eos_3p_hybrid/eos_3p_hybrid.hxx index 7f465b01..d686095f 100644 --- a/EOSX/src/eos_3p_hybrid/eos_3p_hybrid.hxx +++ b/EOSX/src/eos_3p_hybrid/eos_3p_hybrid.hxx @@ -120,6 +120,20 @@ public: return CCTK_REAL(0.0); } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + temp_from_rho_press_ye(const CCTK_REAL rho, CCTK_REAL &press, + const CCTK_REAL ye) const { + printf("eos_3p_hybrid: temperature not implemented"); + return 0.0; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + rho_from_press_temp_ye(CCTK_REAL &press, const CCTK_REAL temp, + const CCTK_REAL ye) const { + printf("eos_3p_hybrid: temperature not implemented"); + return 0.0; + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void press_derivs_from_rho_eps_ye(CCTK_REAL &press, CCTK_REAL &dpdrho, CCTK_REAL &dpdeps, const CCTK_REAL rho, diff --git a/EOSX/src/eos_3p_idealgas/eos_3p_idealgas.hxx b/EOSX/src/eos_3p_idealgas/eos_3p_idealgas.hxx index 96071d25..41250f55 100644 --- a/EOSX/src/eos_3p_idealgas/eos_3p_idealgas.hxx +++ b/EOSX/src/eos_3p_idealgas/eos_3p_idealgas.hxx @@ -76,6 +76,25 @@ public: return temp_over_eps * eps; } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + temp_from_rho_press_ye( + const CCTK_REAL rho, ///< Rest mass density \f$ \rho \f$ + CCTK_REAL &press, ///< Pressure \f$ P \f$ + const CCTK_REAL ye ///< Electron fraction \f$ Y_e \f$ + ) const { + CCTK_REAL eps = eps_from_rho_press_ye(rho, press, ye); + return temp_over_eps * eps; + } + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + rho_from_press_temp_ye( + CCTK_REAL &press, ///< Pressure \f$ P \f$ + const CCTK_REAL temp, ///< Temperature \f$ T \f$ + const CCTK_REAL ye ///< Electron fraction \f$ Y_e \f$ + ) const{ + return press * temp_over_eps / (temp * gm1); + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void press_derivs_from_rho_eps_ye( CCTK_REAL &press, ///< Pressure \f$ P \f$ diff --git a/EOSX/src/eos_3p_tabulated3d/eos_3p_tabulated3d.hxx b/EOSX/src/eos_3p_tabulated3d/eos_3p_tabulated3d.hxx index befeb03e..14c3dbf0 100644 --- a/EOSX/src/eos_3p_tabulated3d/eos_3p_tabulated3d.hxx +++ b/EOSX/src/eos_3p_tabulated3d/eos_3p_tabulated3d.hxx @@ -198,6 +198,37 @@ public: func); } + template + CCTK_HOST CCTK_DEVICE inline CCTK_REAL + logrho_from_var_temp_ye(CCTK_REAL &invar, const CCTK_REAL temp, + const CCTK_REAL ye) const { + // bound inputs + CCTK_REAL t = std::fmin(std::fmax(temp, rgtemp.min), rgtemp.max); + CCTK_REAL lt = std::log(t); + + // table‐edge clamp + auto vmin = + interptable->interpolate(interptable->xmin<0>(), lt, ye)[0]; + auto vmax = + interptable->interpolate(interptable->xmax<0>(), lt, ye)[0]; + if (invar <= vmin) { + invar = vmin; + return interptable->xmin<0>(); + } + if (invar >= vmax) { + invar = vmax; + return interptable->xmax<0>(); + } + + // root‐find for logrho + auto func = [&](CCTK_REAL &lrho) { + CCTK_REAL val = interptable->interpolate(lrho, lt, ye)[0]; + return invar - val; + }; + return zero_brent(interptable->xmin<0>(), interptable->xmax<0>(), 1.e-14, + func); + } + CCTK_HOST CCTK_DEVICE inline CCTK_REAL logtemp_from_rho_eps_ye(const CCTK_REAL rho, CCTK_REAL &eps, const CCTK_REAL ye) const { @@ -223,6 +254,24 @@ public: return exp(lt); } + CCTK_HOST CCTK_DEVICE inline CCTK_REAL + temp_from_rho_press_ye(const CCTK_REAL rho, CCTK_REAL &press, + const CCTK_REAL ye) const { + CCTK_REAL lP = log(press); + CCTK_REAL lt = logtemp_from_rho_var_ye(rho, lP, ye); + press = exp(lP); + return exp(lt); + } + + CCTK_HOST CCTK_DEVICE inline CCTK_REAL + rho_from_press_temp_ye(CCTK_REAL &press, const CCTK_REAL temp, + const CCTK_REAL ye) const { + CCTK_REAL lP = log(press); + CCTK_REAL lr = logrho_from_var_temp_ye(lP, temp, ye); + press = exp(lP); + return exp(lr); + } + CCTK_HOST CCTK_DEVICE inline CCTK_REAL press_from_rho_temp_ye(const CCTK_REAL rho, const CCTK_REAL temp, const CCTK_REAL ye) const {