diff --git a/AsterSeeds/param.ccl b/AsterSeeds/param.ccl index c5907c53..b0cacbd8 100644 --- a/AsterSeeds/param.ccl +++ b/AsterSeeds/param.ccl @@ -2,6 +2,7 @@ KEYWORD test_type "Type of test to set up" STEERABLE=never { + "None" :: "Initial data should be written by another thorn" "1DTest" :: "" "2DTest" :: "" "3DTest" :: "" @@ -149,7 +150,14 @@ REAL dipole_z[2] "z-coordinate of the dipole center" STEERABLE=ALWAYS *:* :: "Anything" } 0.0 +# entropy + +BOOLEAN set_entropy_postinitial "Set the entropy consistently in HydroBaseX_PostInitial" STEERABLE=always +{ +} "yes" + SHARES: EOSX +USES KEYWORD evolution_eos USES CCTK_REAL poly_gamma USES CCTK_REAL poly_k USES CCTK_REAL gl_gamma diff --git a/AsterSeeds/schedule.ccl b/AsterSeeds/schedule.ccl index f07119ec..6bd53885 100644 --- a/AsterSeeds/schedule.ccl +++ b/AsterSeeds/schedule.ccl @@ -91,3 +91,14 @@ if (CCTK_Equals(test_type, "3DTest")) { } } + +#Initial conditions for entropy + +schedule SetEntropy IN HydroBaseX_PostInitial +{ + LANG: C + + READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) + WRITES: HydroBaseX::entropy(everywhere) + +} "Set initial entropy" diff --git a/AsterSeeds/src/SetEntropy.cxx b/AsterSeeds/src/SetEntropy.cxx new file mode 100644 index 00000000..8cf1ac0d --- /dev/null +++ b/AsterSeeds/src/SetEntropy.cxx @@ -0,0 +1,40 @@ +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "eos.hxx" +#include "eos_idealgas.hxx" + +extern "C" void SetEntropy(CCTK_ARGUMENTS) +{ + + using namespace EOSX; + + DECLARE_CCTK_ARGUMENTSX_SetEntropy; + DECLARE_CCTK_PARAMETERS; + + if (set_entropy_postinitial) { + + eos::range rgeps(eps_min, eps_max), rgrho(rho_min, rho_max), + rgye(ye_min, ye_max); + + const eos_idealgas eos_th(gl_gamma, particle_mass, rgeps, rgrho, rgye); + + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, + [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + + entropy(p.I) = eos_th.kappa_from_valid_rho_eps_ye(rho(p.I),eps(p.I),1.0); + + }); + } +} + + diff --git a/AsterSeeds/src/make.code.defn b/AsterSeeds/src/make.code.defn index 864c3f03..cbcfaf32 100644 --- a/AsterSeeds/src/make.code.defn +++ b/AsterSeeds/src/make.code.defn @@ -7,8 +7,8 @@ SRCS = \ 3D_tests.cxx \ atmosphere.cxx \ magTOV.cxx \ - magBNS.cxx - + magBNS.cxx \ + SetEntropy.cxx # Subdirectories containing source files SUBDIRS = diff --git a/AsterX/configuration.ccl b/AsterX/configuration.ccl index 8320a5e3..94b46419 100644 --- a/AsterX/configuration.ccl +++ b/AsterX/configuration.ccl @@ -1,3 +1,7 @@ # Configuration definitions for thorn AsterX REQUIRES Loop EOSX Con2PrimFactory ReconX AsterUtils + +PROVIDES AsterX +{ +} diff --git a/AsterX/interface.ccl b/AsterX/interface.ccl index 3d0e0d49..c6b6c40a 100644 --- a/AsterX/interface.ccl +++ b/AsterX/interface.ccl @@ -8,7 +8,7 @@ USES INCLUDE HEADER: fixmath.hxx USES INCLUDE HEADER: loop_device.hxx USES INCLUDE HEADER: eos_1p.hxx eos_polytropic.hxx USES INCLUDE HEADER: eos.hxx eos_idealgas.hxx -USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx +USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx USES INCLUDE HEADER: reconstruct.hxx USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx @@ -23,11 +23,24 @@ CCTK_REAL mom TYPE=gf CENTERING={ccc} TAGS='rhs="momrhs" dependents="HydroBaseX: CCTK_REAL tau TYPE=gf CENTERING={ccc} TAGS='rhs="taurhs" dependents="HydroBaseX::eps HydroBaseX::press HydroBaseX::temperature HydroBaseX::Ye HydroBaseX::Bvec TmunuBaseX::eTij"' "Conserved internal energy density" +# DEnt is the conserved "evolved" entropy. +# The corresponding primitive is stored in HydroBaseX::entropy. +# However, note that this gf stores the "evolved" entropy which is +# not necessarily the "physical" entropy. This depends on the EOS, e.g. +# for the ideal gas we have entropy = p rho^(-gamma). +# The distinction between "evolved" and "physical" entropy is made +# explicit in EOSX where functions with "entropy_..." and "kappa_..." +# refer to "physical" and "evolved" entropy, respectively. In this +# evolution thorn we always refer to DEnt and entropy to describe +# the "evolved" quantities. +CCTK_REAL DEnt TYPE=gf CENTERING={ccc} TAGS='rhs="DEntrhs" dependents="HydroBaseX::entropy AsterX::flux_x AsterX::flux_y AsterX::flux_z"' "Advected entropy density" + CCTK_REAL flux_x TYPE=gf CENTERING={vcc} TAGS='checkpoint="no"' { fxdens fxmomx fxmomy fxmomz fxtau + fxDEnt fxBx fxBy fxBz } "Fluxes in x direction" @@ -36,6 +49,7 @@ CCTK_REAL flux_y TYPE=gf CENTERING={cvc} TAGS='checkpoint="no"' fydens fymomx fymomy fymomz fytau + fyDEnt fyBx fyBy fyBz } "Fluxes in y direction" @@ -44,6 +58,7 @@ CCTK_REAL flux_z TYPE=gf CENTERING={ccv} TAGS='checkpoint="no"' fzdens fzmomx fzmomy fzmomz fztau + fzDEnt fzBx fzBy fzBz } "Fluxes in z direction" @@ -56,6 +71,7 @@ CCTK_REAL momrhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' CCTK_REAL taurhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Conserved internal energy density RHS" +CCTK_REAL DEntrhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Advected entropy density RHS" CCTK_REAL Avec_x TYPE=gf CENTERING={cvv} TAGS='rhs="Avec_x_rhs"' "x-component of vector potential" CCTK_REAL Avec_y TYPE=gf CENTERING={vcv} TAGS='rhs="Avec_y_rhs"' "y-component of vector potential" diff --git a/AsterX/par/Scaling_magTOV_Z4c_unigrid.par b/AsterX/par/Scaling_magTOV_Z4c_unigrid.par index 160e87e9..d3c06b68 100644 --- a/AsterX/par/Scaling_magTOV_Z4c_unigrid.par +++ b/AsterX/par/Scaling_magTOV_Z4c_unigrid.par @@ -110,7 +110,6 @@ Con2PrimFactory::unit_test = "yes" Con2PrimFactory::B_lim = 1e8 Con2PrimFactory::vw_lim = 1e8 Con2PrimFactory::Ye_lenient = "yes" -Con2PrimFactory::rho_strict = 6.4e-05 EOSX::evolution_eos = "IdealGas" EOSX::gl_gamma = 2.0 diff --git a/AsterX/par/magTOV_Cowling_unigrid.par b/AsterX/par/magTOV_Cowling_unigrid.par index 18d8d9fc..99fa72c8 100644 --- a/AsterX/par/magTOV_Cowling_unigrid.par +++ b/AsterX/par/magTOV_Cowling_unigrid.par @@ -110,7 +110,6 @@ Con2PrimFactory::unit_test = "yes" Con2PrimFactory::B_lim = 1e8 Con2PrimFactory::vw_lim = 1e8 Con2PrimFactory::Ye_lenient = "yes" -Con2PrimFactory::rho_strict = 6.4e-05 EOSX::evolution_eos = "IdealGas" EOSX::gl_gamma = 2.0 diff --git a/AsterX/par/magTOV_Z4c_AMR.par b/AsterX/par/magTOV_Z4c_AMR.par index 9cb9bb28..fabf226a 100644 --- a/AsterX/par/magTOV_Z4c_AMR.par +++ b/AsterX/par/magTOV_Z4c_AMR.par @@ -121,7 +121,6 @@ Con2PrimFactory::unit_test = "yes" Con2PrimFactory::B_lim = 1e8 Con2PrimFactory::vw_lim = 1e8 Con2PrimFactory::Ye_lenient = "yes" -Con2PrimFactory::rho_strict = 6.4e-05 EOSX::evolution_eos = "IdealGas" EOSX::gl_gamma = 2.0 diff --git a/AsterX/param.ccl b/AsterX/param.ccl index b3cad544..71c46f38 100644 --- a/AsterX/param.ccl +++ b/AsterX/param.ccl @@ -15,6 +15,10 @@ BOOLEAN use_uct "Shall we use the Upwind-CT method to compute the electric field { } "no" +BOOLEAN use_entropy_fix "Shall we use inversion based on entropy as a backup?" STEERABLE=always +{ +} "no" + KEYWORD flux_type "Flux solver" STEERABLE=always { "LxF" :: "" @@ -108,7 +112,6 @@ USES CCTK_REAL p_atmo USES CCTK_REAL Ye_atmo USES CCTK_REAL B_lim USES CCTK_REAL vw_lim -USES CCTK_REAL rho_strict USES BOOLEAN Ye_lenient USES CCTK_INT max_iter USES CCTK_REAL c2p_tol @@ -120,6 +123,8 @@ USES CCTK_REAL alp_thresh USES CCTK_REAL rho_BH USES CCTK_REAL eps_BH USES CCTK_REAL vwlim_BH +USES CCTK_REAL cons_error_limit +USES BOOLEAN use_z SHARES: ReconX USES KEYWORD reconstruction_method diff --git a/AsterX/schedule.ccl b/AsterX/schedule.ccl index 28db076e..8b67a204 100644 --- a/AsterX/schedule.ccl +++ b/AsterX/schedule.ccl @@ -1,8 +1,8 @@ # Schedule definitions for thorn AsterX -STORAGE: dens mom tau dB Psi HydroBaseX::Bvec dBx_stag dBy_stag dBz_stag +STORAGE: dens mom tau DEnt dB Psi HydroBaseX::Bvec dBx_stag dBy_stag dBz_stag STORAGE: flux_x flux_y flux_z -STORAGE: densrhs momrhs taurhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs +STORAGE: densrhs momrhs taurhs DEntrhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs STORAGE: ADMBaseX::metric ADMBaseX::lapse ADMBaseX::shift ADMBaseX::curv STORAGE: Aux_in_RHSof_A_Psi STORAGE: TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij @@ -69,12 +69,13 @@ SCHEDULE AsterX_Prim2Con_Initial IN AsterX_InitialGroup AFTER AsterX_ComputeBFro LANG: C READS: ADMBaseX::metric(interior) READS: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior) - WRITES: dens(interior) tau(interior) mom(interior) dB(interior) + READS: HydroBaseX::entropy(interior) + WRITES: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior) WRITES: Psi(everywhere) WRITES: saved_prims WRITES: zvec WRITES: svec - SYNC: dens tau mom dB + SYNC: dens tau DEnt mom dB SYNC: saved_prims SYNC: zvec SYNC: svec @@ -86,7 +87,7 @@ SCHEDULE AsterX_Sync AT postregrid { LANG: C OPTIONS: global - SYNC: dens tau mom Avec_x Avec_y Avec_z Psi + SYNC: dens tau mom DEnt Avec_x Avec_y Avec_z Psi SYNC: saved_prims } "Synchronize" @@ -94,7 +95,7 @@ SCHEDULE AsterX_Sync IN ODESolvers_PostStep { LANG: C OPTIONS: global - SYNC: dens tau mom Avec_x Avec_y Avec_z Psi + SYNC: dens tau mom DEnt Avec_x Avec_y Avec_z Psi } "Synchronize" @@ -124,21 +125,23 @@ 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: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior) READS: saved_prims(interior) READS: Avec_x(interior) Avec_y(interior) Avec_z(interior) WRITES: con2prim_flag(interior) - WRITES: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior) + WRITES: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior) + WRITES: HydroBaseX::entropy(interior) WRITES: saved_prims(interior) WRITES: zvec(interior) WRITES: svec(interior) - WRITES: dens(interior) tau(interior) mom(interior) dB(interior) + WRITES: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior) SYNC: con2prim_flag - SYNC: HydroBaseX::rho HydroBaseX::vel HydroBaseX::eps HydroBaseX::press HydroBaseX::Bvec + SYNC: HydroBaseX::rho HydroBaseX::vel HydroBaseX::eps HydroBaseX::press HydroBaseX::Bvec + SYNC: HydroBaseX::entropy SYNC: saved_prims SYNC: zvec SYNC: svec - SYNC: dens tau mom dB + SYNC: dens tau DEnt mom dB } "Calculate primitive variables from conservative variables" @@ -153,9 +156,10 @@ SCHEDULE AsterX_Fluxes IN AsterX_RHSGroup READS: ADMBaseX::metric(everywhere) READS: ADMBaseX::lapse(everywhere) READS: ADMBaseX::shift(everywhere) - READS: dens(everywhere) tau(everywhere) mom(everywhere) + READS: dens(everywhere) tau(everywhere) mom(everywhere) DEnt(everywhere) READS: HydroBaseX::rho(everywhere) HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) READS: HydroBaseX::Bvec(everywhere) + READS: HydroBaseX::entropy(everywhere) READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere) READS: zvec_x(everywhere) zvec_y(everywhere) zvec_z(everywhere) READS: svec_x(everywhere) svec_y(everywhere) svec_z(everywhere) @@ -181,8 +185,8 @@ 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: densrhs(interior) taurhs(interior) momrhs(interior) - SYNC: densrhs taurhs momrhs + WRITES: densrhs(interior) taurhs(interior) momrhs(interior) DEntrhs(interior) + SYNC: densrhs taurhs momrhs DEntrhs } "Calculate the source terms and compute the RHS of the hydro equations" SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms @@ -191,15 +195,15 @@ SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere) READS: HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere) - READS: densrhs(everywhere) taurhs(everywhere) momrhs(everywhere) + READS: densrhs(everywhere) taurhs(everywhere) momrhs(everywhere) DEntrhs(everywhere) READS: Psi(everywhere) READS: Aux_in_RHSof_A_Psi(everywhere) READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere) READS: vtilde_xface(everywhere) vtilde_yface(everywhere) vtilde_zface(everywhere) READS: a_xface(everywhere) a_yface(everywhere) a_zface(everywhere) - WRITES: densrhs(interior) taurhs(interior) momrhs(interior) + WRITES: densrhs(interior) taurhs(interior) momrhs(interior) DEntrhs(interior) WRITES: Avec_x_rhs(interior) Avec_y_rhs(interior) Avec_z_rhs(interior) Psi_rhs(interior) - SYNC: densrhs taurhs momrhs + SYNC: densrhs taurhs momrhs DEntrhs SYNC: Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs } "Update the RHS of the hydro equations with the flux contributions" diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 0d4ca3e1..5fc2d4fd 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -8,6 +8,7 @@ #include "c2p.hxx" #include "c2p_1DPalenzuela.hxx" #include "c2p_2DNoble.hxx" +#include "c2p_1DEntropy.hxx" #include "eos_1p.hxx" #include "eos_polytropic.hxx" @@ -59,6 +60,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, 1, 1, 1>(grid.nghostzones, [=] CCTK_DEVICE( const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // Note that HydroBaseX gfs are NaN when entering this loop due + // explicit dependence on conservatives from + // AsterX -> dependents tag + // Setting up atmosphere CCTK_REAL rho_atm = 0.0; // dummy initialization CCTK_REAL press_atm = 0.0; // dummy initialization @@ -83,15 +88,26 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, 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); + CCTK_REAL entropy_atm = eos_th.kappa_from_valid_rho_eps_ye(rho_atm, eps_atm, Ye_atmo); + atmosphere atmo(rho_atm, eps_atm, Ye_atmo, press_atm, entropy_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); + c2p_2DNoble c2p_Noble(eos_th, atmo, max_iter, c2p_tol, + alp_thresh, cons_error_limit, + vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, + Ye_lenient, use_z); // Construct Palenzuela c2p object: - c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, rho_strict, - vw_lim, B_lim, Ye_lenient); + c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, + alp_thresh, cons_error_limit, + vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, + Ye_lenient, use_z); + + // Construct Entropy c2p object: + c2p_1DEntropy c2p_Ent(eos_th, atmo, max_iter, c2p_tol, + alp_thresh, cons_error_limit, + vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, + Ye_lenient, use_z); /* Get covariant metric */ const smat glo( @@ -102,8 +118,24 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, const CCTK_REAL sqrt_detg = sqrt(spatial_detg); vec v_up{saved_velx(p.I), saved_vely(p.I), saved_velz(p.I)}; - const vec v_low = calc_contraction(glo, v_up); - CCTK_REAL wlor = calc_wlorentz(v_low, v_up); + vec v_low = calc_contraction(glo, v_up); + CCTK_REAL zsq{0.0}; + + // TODO: Debug code to capture v>1 early, + // remove soon + const CCTK_REAL vsq = calc_contraction(v_low,v_up); + if (vsq >= 1.0) { + CCTK_REAL wlim = sqrt(1.0 + vw_lim * vw_lim); + CCTK_REAL vlim = vw_lim/wlim; + v_up *= vlim/sqrt(vsq); + v_low *= vlim/sqrt(vsq); + zsq = vw_lim; + } else { + zsq = vsq/(1.0-vsq); + } + + //CCTK_REAL wlor = calc_wlorentz(v_low, v_up); + CCTK_REAL wlor = sqrt(1.0+zsq); vec Bup{dBx(p.I) / sqrt_detg, dBy(p.I) / sqrt_detg, dBz(p.I) / sqrt_detg}; @@ -111,13 +143,15 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, CCTK_REAL dummy_Ye = 0.5; CCTK_REAL dummy_dYe = 0.5; prim_vars pv; - prim_vars pv_seeds{saved_rho(p.I), saved_eps(p.I), dummy_Ye, press(p.I), + prim_vars pv_seeds{saved_rho(p.I), saved_eps(p.I), dummy_Ye, press(p.I), entropy(p.I), v_up, wlor, Bup}; + // Note that pv_seeds.press and pv_seeds.entropy are NaN at this point // Note that cv are densitized, i.e. they all include sqrt_detg cons_vars cv{dens(p.I), {momx(p.I), momy(p.I), momz(p.I)}, tau(p.I), dummy_dYe, + DEnt(p.I), {dBx(p.I), dBy(p.I), dBz(p.I)}}; if (dens(p.I) <= sqrt_detg * rho_atmo_cut) { @@ -136,28 +170,17 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, 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); - // check on velocities - CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); - CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; - CCTK_REAL sol_v = - sqrt((pv_seeds.w_lor * pv_seeds.w_lor - 1.0)) / pv_seeds.w_lor; - if (sol_v > vlim_BH) { - pv_seeds.vel *= vlim_BH / sol_v; - pv_seeds.w_lor = wlim_BH; - } - cv.from_prim(pv_seeds, glo); + c2p_Noble.bh_interior_fail(eos_th,pv,cv,glo); } } // Construct error report object: c2p_report rep_first; c2p_report rep_second; + c2p_report rep_ent; + + /* set flag to success */ + con2prim_flag(p.I) = 1; // Calling the first C2P switch (c2p_fir) { @@ -166,7 +189,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, break; } case c2p_first_t::Palenzuela: { - c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, rep_first); + //c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, rep_first); + c2p_Pal.solve(eos_th, pv, cv, glo, rep_first); break; } default: @@ -174,10 +198,10 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, } if (rep_first.failed()) { - if (debug_mode) { - printf("First C2P failed :( \n"); - rep_first.debug_message(); - printf("Calling the back up C2P.. \n"); + if (debug_mode) { + printf("First C2P failed :( \n"); + rep_first.debug_message(); + printf("Calling the back up C2P.. \n"); } // Calling the second C2P switch (c2p_sec) { @@ -186,7 +210,8 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, break; } case c2p_second_t::Palenzuela: { - c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, rep_second); + //c2p_Pal.solve(eos_th, pv, pv_seeds, cv, glo, rep_second); + c2p_Pal.solve(eos_th, pv, cv, glo, rep_second); break; } default: @@ -195,83 +220,103 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, } if (rep_first.failed() && rep_second.failed()) { - if (debug_mode) { - 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); - // check on velocities - CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); - CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; - CCTK_REAL sol_v = sqrt((pv.w_lor * pv.w_lor - 1.0)) / pv.w_lor; - if (sol_v > vlim_BH) { - pv.vel *= vlim_BH / sol_v; - pv.w_lor = wlim_BH; + + if (use_entropy_fix) { + + c2p_Ent.solve(eos_th, pv, cv, glo, rep_ent); + + if (rep_ent.failed()) { + + con2prim_flag(p.I) = 0; + + if (debug_mode) { + printf("Entropy C2P failed. Setting point to atmosphere.\n"); + rep_ent.debug_message(); + printf( + "WARNING: \n" + "C2Ps failed. Printing cons and saved prims before set to " + "atmo: \n" + "cctk_iteration = %i \n " + "x, y, z = %26.16e, %26.16e, %26.16e \n " + "dens = %26.16e \n tau = %26.16e \n momx = %26.16e \n " + "momy = %26.16e \n momz = %26.16e \n dBx = %26.16e \n " + "dBy = %26.16e \n dBz = %26.16e \n " + "saved_rho = %26.16e \n saved_eps = %26.16e \n press= %26.16e \n " + "saved_velx = %26.16e \n saved_vely = %26.16e \n saved_velz = " + "%26.16e \n " + "Bvecx = %26.16e \n Bvecy = %26.16e \n " + "Bvecz = %26.16e \n " + "Avec_x = %26.16e \n Avec_y = %26.16e \n Avec_z = %26.16e \n ", + cctk_iteration, p.x, p.y, p.z, dens(p.I), tau(p.I), momx(p.I), + momy(p.I), momz(p.I), dBx(p.I), dBy(p.I), dBz(p.I), pv.rho, pv.eps, + pv.press, pv.vel(0), pv.vel(1), pv.vel(2), pv.Bvec(0), pv.Bvec(1), + pv.Bvec(2), + // rho(p.I), eps(p.I), press(p.I), velx(p.I), vely(p.I), + // velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), + Avec_x(p.I), Avec_y(p.I), Avec_z(p.I)); } - cv.from_prim(pv, glo); - rep_first.set_atmo = 0; - rep_second.set_atmo = 0; - } - } - con2prim_flag(p.I) = 0; - } - if (rep_first.set_atmo && rep_second.set_atmo) { - if (debug_mode) { - printf( - "WARNING: \n" - "C2Ps failed. Printing cons and saved prims before set to " - "atmo: \n" - "cctk_iteration = %i \n " - "x, y, z = %26.16e, %26.16e, %26.16e \n " - "dens = %26.16e \n tau = %26.16e \n momx = %26.16e \n " - "momy = %26.16e \n momz = %26.16e \n dBx = %26.16e \n " - "dBy = %26.16e \n dBz = %26.16e \n " - "saved_rho = %26.16e \n saved_eps = %26.16e \n press= %26.16e \n " - "saved_velx = %26.16e \n saved_vely = %26.16e \n saved_velz = " - "%26.16e \n " - "Bvecx = %26.16e \n Bvecy = %26.16e \n " - "Bvecz = %26.16e \n " - "Avec_x = %26.16e \n Avec_y = %26.16e \n Avec_z = %26.16e \n ", - cctk_iteration, p.x, p.y, p.z, dens(p.I), tau(p.I), momx(p.I), - momy(p.I), momz(p.I), dBx(p.I), dBy(p.I), dBz(p.I), pv.rho, pv.eps, - pv.press, pv.vel(0), pv.vel(1), pv.vel(2), pv.Bvec(0), pv.Bvec(1), - pv.Bvec(2), - // rho(p.I), eps(p.I), press(p.I), velx(p.I), vely(p.I), - // velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), - Avec_x(p.I), Avec_y(p.I), Avec_z(p.I)); - } + if ( (alp(p.I) < alp_thresh) ) { + c2p_Noble.bh_interior_fail(eos_th,pv,cv,glo); + } else { + // set to atmo + cv.dBvec(0) = dBx(p.I); + cv.dBvec(1) = dBy(p.I); + cv.dBvec(2) = dBz(p.I); + pv.Bvec = cv.dBvec / sqrt_detg; + atmo.set(pv, cv, glo); + } + } - // set to atmo - cv.dBvec(0) = dBx(p.I); - cv.dBvec(1) = dBy(p.I); - cv.dBvec(2) = dBz(p.I); - pv.Bvec = cv.dBvec / sqrt_detg; - atmo.set(pv, cv, glo); + } else { + + con2prim_flag(p.I) = 0; + + if (debug_mode) { + printf("Second C2P failed too :( :( \n"); + rep_second.debug_message(); + printf( + "WARNING: \n" + "C2Ps failed. Printing cons and saved prims before set to " + "atmo: \n" + "cctk_iteration = %i \n " + "x, y, z = %26.16e, %26.16e, %26.16e \n " + "dens = %26.16e \n tau = %26.16e \n momx = %26.16e \n " + "momy = %26.16e \n momz = %26.16e \n dBx = %26.16e \n " + "dBy = %26.16e \n dBz = %26.16e \n " + "saved_rho = %26.16e \n saved_eps = %26.16e \n press= %26.16e \n " + "saved_velx = %26.16e \n saved_vely = %26.16e \n saved_velz = " + "%26.16e \n " + "Bvecx = %26.16e \n Bvecy = %26.16e \n " + "Bvecz = %26.16e \n " + "Avec_x = %26.16e \n Avec_y = %26.16e \n Avec_z = %26.16e \n ", + cctk_iteration, p.x, p.y, p.z, dens(p.I), tau(p.I), momx(p.I), + momy(p.I), momz(p.I), dBx(p.I), dBy(p.I), dBz(p.I), pv.rho, pv.eps, + pv.press, pv.vel(0), pv.vel(1), pv.vel(2), pv.Bvec(0), pv.Bvec(1), + pv.Bvec(2), + // rho(p.I), eps(p.I), press(p.I), velx(p.I), vely(p.I), + // velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), + Avec_x(p.I), Avec_y(p.I), Avec_z(p.I)); + } - // assert(0); + if ( (alp(p.I) < alp_thresh) && ( (pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH) ) ) { + c2p_Noble.bh_interior_fail(eos_th,pv,cv,glo); + } else { + // set to atmo + cv.dBvec(0) = dBx(p.I); + cv.dBvec(1) = dBy(p.I); + cv.dBvec(2) = dBz(p.I); + pv.Bvec = cv.dBvec / sqrt_detg; + atmo.set(pv, cv, glo); + } + } } - /* set flag to success */ - con2prim_flag(p.I) = 1; - // dummy vars CCTK_REAL Ex, Ey, Ez; // Write back pv - pv.scatter(rho(p.I), eps(p.I), dummy_Ye, press(p.I), velx(p.I), vely(p.I), + pv.scatter(rho(p.I), eps(p.I), dummy_Ye, press(p.I), entropy(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); @@ -287,7 +332,7 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, // Write back cv cv.scatter(dens(p.I), momx(p.I), momy(p.I), momz(p.I), tau(p.I), dummy_Ye, - dBx(p.I), dBy(p.I), dBz(p.I)); + DEnt(p.I), dBx(p.I), dBy(p.I), dBz(p.I)); // Update saved prims saved_rho(p.I) = rho(p.I); @@ -295,6 +340,7 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, saved_vely(p.I) = vely(p.I); saved_velz(p.I) = velz(p.I); saved_eps(p.I) = eps(p.I); + }); // Loop } diff --git a/AsterX/src/con2prim_rpa.cxx b/AsterX/src/con2prim_rpa.cxx index 5fc50b0c..d9b5dff8 100644 --- a/AsterX/src/con2prim_rpa.cxx +++ b/AsterX/src/con2prim_rpa.cxx @@ -74,7 +74,7 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { CCTK_REAL dummy_dYe = 0.5; // Get a recovery function - con2prim_mhd cv2pv(eos, rho_strict, Ye_lenient, vw_lim, B_lim, atmo, + con2prim_mhd cv2pv(eos, 1.0e5, Ye_lenient, vw_lim, B_lim, atmo, c2p_tol, max_iter); /* Get covariant metric */ @@ -109,6 +109,9 @@ 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)}}; + // Invert entropy here + entropy(p.I) = DEnt(p.I)/dens(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 diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx index 1fe0ba84..7f14092c 100644 --- a/AsterX/src/fluxes.cxx +++ b/AsterX/src/fluxes.cxx @@ -44,6 +44,7 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { /* grid functions for fluxes */ const vec, dim> fluxdenss{fxdens, fydens, fzdens}; + const vec, dim> fluxDEnts{fxDEnt, fyDEnt, fzDEnt}; const vec, dim> fluxmomxs{fxmomx, fymomx, fzmomx}; const vec, dim> fluxmomys{fxmomy, fymomy, fzmomy}; const vec, dim> fluxmomzs{fxmomz, fymomz, fzmomz}; @@ -214,6 +215,8 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { vec rho_rc{reconstruct_pt(rho, p, true, true)}; + vec entropy_rc{reconstruct_pt(entropy, p, true, true)}; + // set to atmo if reconstructed rho is less than atmo or is negative if (rho_rc(0) < rho_abs_min) { rho_rc(0) = rho_abs_min; @@ -387,6 +390,12 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { return sqrtg * rho_rc(f) * w_lorentz_rc(f); }); + /* DEnt = sqrt(g) * D * s = sqrt(g) * (rho * W) * s */ + /* s = entropy */ + const vec DEnt_rc([&](int f) ARITH_INLINE { + return sqrtg * rho_rc(f) * w_lorentz_rc(f) * entropy_rc(f); + }); + /* auxiliary: dens * h * W = sqrt(g) * rho * h * W^2 */ const vec dens_h_W_rc([&](int f) ARITH_INLINE { return dens_rc(f) * h_rc(f) * w_lorentz_rc(f); @@ -428,6 +437,10 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { const vec flux_dens( [&](int f) ARITH_INLINE { return dens_rc(f) * vtilde_rc(f); }); + /* flux(DEnt) = sqrt(g) * D * s * vtilde^i = sqrt(g) * rho * W * s * vtilde^i */ + const vec flux_DEnt( + [&](int f) ARITH_INLINE { return DEnt_rc(f) * vtilde_rc(f); }); + /* flux(mom_j)^i = sqrt(g)*( * S_j*vtilde^i + alpha*((pgas+pmag)*delta^i_j - b_jB^i/W) ) */ const vec, 3> flux_moms([&](int j) ARITH_INLINE { @@ -464,6 +477,7 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) { /* Calculate numerical fluxes */ fluxdenss(dir)(p.I) = calcflux(lambda, dens_rc, flux_dens); + fluxDEnts(dir)(p.I) = calcflux(lambda, DEnt_rc, flux_DEnt); fluxmomxs(dir)(p.I) = calcflux(lambda, moms_rc(0), flux_moms(0)); fluxmomys(dir)(p.I) = calcflux(lambda, moms_rc(1), flux_moms(1)); fluxmomzs(dir)(p.I) = calcflux(lambda, moms_rc(2), flux_moms(2)); diff --git a/AsterX/src/prim2con.cxx b/AsterX/src/prim2con.cxx index 0f435d76..27092968 100644 --- a/AsterX/src/prim2con.cxx +++ b/AsterX/src/prim2con.cxx @@ -44,6 +44,8 @@ extern "C" void AsterX_Prim2Con_Initial(CCTK_ARGUMENTS) { dBy(p.I) = cv.dBvec(1); dBz(p.I) = cv.dBvec(2); + DEnt(p.I) = entropy(p.I)*cv.dens; + saved_rho(p.I) = pv.rho; saved_velx(p.I) = pv.vel(0); saved_vely(p.I) = pv.vel(1); diff --git a/AsterX/src/rhs.cxx b/AsterX/src/rhs.cxx index ba161b32..562b684d 100644 --- a/AsterX/src/rhs.cxx +++ b/AsterX/src/rhs.cxx @@ -71,6 +71,7 @@ extern "C" void AsterX_RHS(CCTK_ARGUMENTS) { 1 / CCTK_DELTA_SPACE(2)}; const vec, dim> gf_fdens{fxdens, fydens, fzdens}; + const vec, dim> gf_fDEnt{fxDEnt, fyDEnt, fzDEnt}; const vec, dim> gf_fmomx{fxmomx, fymomx, fzmomx}; const vec, dim> gf_fmomy{fxmomy, fymomy, fzmomy}; const vec, dim> gf_fmomz{fxmomz, fymomz, fzmomz}; @@ -187,6 +188,7 @@ extern "C" void AsterX_RHS(CCTK_ARGUMENTS) { grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { densrhs(p.I) += calcupdate_hydro(gf_fdens, p); + DEntrhs(p.I) += calcupdate_hydro(gf_fDEnt, p); momxrhs(p.I) += calcupdate_hydro(gf_fmomx, p); momyrhs(p.I) += calcupdate_hydro(gf_fmomy, p); momzrhs(p.I) += calcupdate_hydro(gf_fmomz, p); diff --git a/AsterX/src/source.cxx b/AsterX/src/source.cxx index d27487c5..7ded0a5f 100644 --- a/AsterX/src/source.cxx +++ b/AsterX/src/source.cxx @@ -159,6 +159,7 @@ template void SourceTerms(CCTK_ARGUMENTS) { /* Update the RHS grid functions */ densrhs(p.I) = 0.0; + DEntrhs(p.I) = 0.0; momxrhs(p.I) = alp_avg * sqrt_detg * mom_source(0); momyrhs(p.I) = alp_avg * sqrt_detg * mom_source(1); momzrhs(p.I) = alp_avg * sqrt_detg * mom_source(2); diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl index 0470f717..105da7fd 100644 --- a/Con2PrimFactory/interface.ccl +++ b/Con2PrimFactory/interface.ccl @@ -5,7 +5,8 @@ INCLUDES HEADER: atmo.hxx IN atmo.hxx INCLUDES HEADER: prims.hxx IN prims.hxx INCLUDES HEADER: cons.hxx IN cons.hxx INCLUDES HEADER: c2p_2DNoble.hxx IN c2p_2DNoble.hxx -INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx +INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx +INCLUDES HEADER: c2p_1DEntropy.hxx IN c2p_1DEntropy.hxx USES INCLUDE HEADER: eos.hxx USES INCLUDE HEADER: eos_idealgas.hxx USES INCLUDE HEADER: roots.hxx diff --git a/Con2PrimFactory/param.ccl b/Con2PrimFactory/param.ccl index 36b2ec39..b0d8b179 100644 --- a/Con2PrimFactory/param.ccl +++ b/Con2PrimFactory/param.ccl @@ -52,11 +52,6 @@ CCTK_REAL vw_lim "Upper limit for the value of velocity * lorentz_factor " STEER 0: :: "Must be positive" } 10.0 -CCTK_REAL rho_strict "Density above which RePrimAnd c2p corrections are forbidden" STEERABLE=ALWAYS -{ - 0: :: "Must be positive" -} 1e100 - BOOLEAN Ye_lenient "Whether to allow restricting Y_e to a valid range" STEERABLE=ALWAYS { } "no" @@ -113,3 +108,14 @@ CCTK_REAL vwlim_BH " Maximum zvec (v*w) allowed in regions where alp 1 + +BOOLEAN use_z "" STEERABLE=ALWAYS +{ +} "no" diff --git a/Con2PrimFactory/src/atmo.hxx b/Con2PrimFactory/src/atmo.hxx index 35821e64..36c7fce5 100644 --- a/Con2PrimFactory/src/atmo.hxx +++ b/Con2PrimFactory/src/atmo.hxx @@ -18,15 +18,16 @@ struct atmosphere { CCTK_REAL eps_atmo; CCTK_REAL ye_atmo; CCTK_REAL press_atmo; + CCTK_REAL entropy_atmo; CCTK_REAL rho_cut; CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline atmosphere() = default; CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline atmosphere( - CCTK_REAL rho_, CCTK_REAL eps_, CCTK_REAL Ye_, CCTK_REAL press_, + CCTK_REAL rho_, CCTK_REAL eps_, CCTK_REAL Ye_, CCTK_REAL press_, CCTK_REAL entropy_, CCTK_REAL rho_cut_) - : rho_atmo(rho_), eps_atmo(eps_), ye_atmo(Ye_), press_atmo(press_), + : rho_atmo(rho_), eps_atmo(eps_), ye_atmo(Ye_), press_atmo(press_), entropy_atmo(entropy_), rho_cut(rho_cut_) {} CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline atmosphere & @@ -38,6 +39,7 @@ struct atmosphere { eps_atmo = other.eps_atmo; ye_atmo = other.ye_atmo; press_atmo = other.press_atmo; + entropy_atmo = other.entropy_atmo; rho_cut = other.rho_cut; return *this; } @@ -49,6 +51,7 @@ struct atmosphere { pv.eps = eps_atmo; pv.Ye = ye_atmo; pv.press = press_atmo; + pv.entropy = entropy_atmo; pv.vel(0) = 0.0; pv.vel(1) = 0.0; pv.vel(2) = 0.0; @@ -68,6 +71,7 @@ struct atmosphere { cv.mom(1) = 0.0; cv.mom(2) = 0.0; cv.dYe = cv.dens * ye_atmo; + cv.DEnt = cv.dens * entropy_atmo; const vec &B_up = pv.Bvec; const vec B_low = calc_contraction(g, B_up); CCTK_REAL Bsq = calc_contraction(B_up, B_low); diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx index c66e7f67..ec7bdacf 100644 --- a/Con2PrimFactory/src/c2p.hxx +++ b/Con2PrimFactory/src/c2p.hxx @@ -37,29 +37,170 @@ class c2p { protected: /* The constructor must initialize the following variables */ + atmosphere atmo; CCTK_INT maxIterations; CCTK_REAL tolerance; - CCTK_REAL rho_strict; - bool ye_lenient; - CCTK_REAL v_lim; - CCTK_REAL w_lim; + CCTK_REAL alp_thresh; + CCTK_REAL cons_error; CCTK_REAL vw_lim; + CCTK_REAL w_lim; + CCTK_REAL v_lim; CCTK_REAL Bsq_lim; - atmosphere atmo; + CCTK_REAL rho_BH; + CCTK_REAL eps_BH; + CCTK_REAL vwlim_BH; + bool ye_lenient; + bool use_zprim; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Ssq_Exact(vec &mom, + get_Ssq_Exact(const vec &mom, const smat &gup) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Bsq_Exact(vec &B_up, + get_Bsq_Exact(const vec &B_up, const smat &glo) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_BiSi_Exact(vec &Bvec, vec &mom) const; + get_BiSi_Exact(const vec &Bvec, const vec &mom) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec - get_WLorentz_bsq_Seeds(vec &B_up, vec &v_up, + get_WLorentz_bsq_Seeds(const vec &B_up, const vec &v_up, const smat &glo) const; + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + prims_floors_and_ceilings(const EOSType &eos_th, prim_vars &pv, const cons_vars &cv, + const smat &glo, c2p_report &rep) const; + + public: + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + bh_interior_fail(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, + const smat &glo) const; }; +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +c2p::prims_floors_and_ceilings(const EOSType &eos_th, prim_vars &pv, const cons_vars &cv, + const smat &glo, c2p_report &rep) const { + + //const CCTK_REAL spatial_detg = calc_det(glo); + //const CCTK_REAL sqrt_detg = sqrt(spatial_detg); + + // Lower velocity + const vec v_low = calc_contraction(glo, pv.vel); + + // ---------- + // Floor and ceiling for rho and velocity + // Keeps pressure the same and changes eps + // ---------- + + // check if computed velocities are within the specified limit + CCTK_REAL vsq_Sol = calc_contraction(v_low, pv.vel); + CCTK_REAL sol_v = sqrt(vsq_Sol); + if (sol_v > v_lim) { + /* + printf("(sol_v > v_lim) is true! \n"); + printf("sol_v, v_lim: %26.16e, %26.16e \n", sol_v, v_lim); + */ + // add mass, keeps conserved density D + pv.rho = cv.dens / w_lim; + pv.eps = eos_th.eps_from_valid_rho_press_ye(pv.rho, pv.press, pv.Ye); + pv.entropy = + eos_th.kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + // if (pv.rho >= rho_strict) { + // rep.set_speed_limit({ sol_v, sol_v, sol_v }); + // set_to_nan(pv, cv); + // return; + //} + pv.vel *= v_lim / sol_v; + pv.w_lor = w_lim; + // pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), + // eos_th.rgeps.max); + // pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + + rep.adjust_cons = true; + } + + if (pv.rho > eos_th.rgrho.max) { + + // remove mass, changes conserved density D + pv.rho = eos_th.rgrho.max; + pv.eps = eos_th.eps_from_valid_rho_press_ye(eos_th.rgrho.max, pv.press, pv.Ye); + pv.entropy = + eos_th.kappa_from_valid_rho_eps_ye(eos_th.rgrho.max, pv.eps, pv.Ye); + + rep.adjust_cons = true; + } + + // ---------- + // Floor and ceiling for eps + // Keeps rho the same and changes press + // ---------- + + // check the validity of the computed eps + auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); + if (pv.eps > rgeps.max) { + // printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n"); + // if (pv.rho >= rho_strict) { + // rep.set_range_eps(pv.eps); // sets adjust_cons to false by default + // rep.adjust_cons = true; + // set_to_nan(pv, cv); + // return; + //} + pv.eps = rgeps.max; + pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_th.kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + + rep.adjust_cons = true; + } else if (pv.eps < rgeps.min) { + /* + printf( + "(pv.eps < rgeps.min) is true! pv.eps, rgeps.min: %26.16e, %26.16e + \n", + pv.eps, rgeps.min); + printf(" Not adjusting cons.. \n"); + */ + // rep.set_range_eps(rgeps.min); // sets adjust_cons to true + pv.eps = rgeps.min; + pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_th.kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + + rep.adjust_cons = true; + } + + // TODO: check validity for Ye + +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +c2p::bh_interior_fail(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, + const smat &glo) const { + + // 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 + 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 = 0.5; + pv.press = + eos_th.press_from_valid_rho_eps_ye(rho_BH, eps_BH, 0.5); + pv.entropy = + eos_th.kappa_from_valid_rho_eps_ye(rho_BH, eps_BH, 0.5); + // check on velocities + CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); + CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; + CCTK_REAL sol_v = sqrt((pv.w_lor * pv.w_lor - 1.0)) / pv.w_lor; + if (sol_v > vlim_BH) { + pv.vel *= vlim_BH / sol_v; + pv.w_lor = wlim_BH; + } + cv.from_prim(pv, glo); + +} + } // namespace Con2PrimFactory #endif diff --git a/Con2PrimFactory/src/c2p_1DEntropy.hxx b/Con2PrimFactory/src/c2p_1DEntropy.hxx new file mode 100644 index 00000000..fdbaf1fd --- /dev/null +++ b/Con2PrimFactory/src/c2p_1DEntropy.hxx @@ -0,0 +1,441 @@ +#ifndef C2P_1DENTROPY_HXX +#define C2P_1DENTROPY_HXX + +#include "c2p.hxx" +#include "roots.hxx" +namespace Con2PrimFactory { + +class c2p_1DEntropy : public c2p { +public: + /* Some attributes */ + CCTK_REAL GammaIdealFluid; + + /* Constructor */ + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DEntropy( + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z); + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + get_Ssq_Exact(const vec &mom, const smat &gup) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + get_Bsq_Exact(const vec &B_up, const smat &glo) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + get_BiSi_Exact(const vec &Bvec, const vec &mom) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec + get_WLorentz_vsq_bsq_Seeds(const vec &B_up, const vec &v_up, + const smat &glo) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + set_to_nan(prim_vars &pv, cons_vars &cv) const; + /* Called by 1DEntropy */ + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + xEntropyToPrim(CCTK_REAL xEntropy_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq, + CCTK_REAL BiSi, const EOSType &eos_th, prim_vars &pv, const cons_vars &cv, + const smat &gup, + const smat &glo) const; + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + funcRoot_1DEntropy(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, CCTK_REAL x, + const EOSType &eos_th, const cons_vars &cv) const; + + template + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void + solve(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, + const smat &glo, c2p_report &rep) const; + + /* Destructor */ + CCTK_HOST CCTK_DEVICE ~c2p_1DEntropy(); +}; + +/* Constructor */ +template +CCTK_HOST CCTK_DEVICE +CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DEntropy::c2p_1DEntropy( + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z) { + + // Base + atmo = atm; + maxIterations = maxIter; + tolerance = tol; + alp_thresh = alp_thresh_in; + cons_error = consError; + vw_lim = vwlim; + w_lim = sqrt(1.0 + vw_lim * vw_lim); + v_lim = vw_lim / w_lim; + Bsq_lim = B_lim * B_lim; + rho_BH = rho_BH_in; + eps_BH = eps_BH_in; + vwlim_BH = vwlim_BH_in; + ye_lenient = ye_len; + use_zprim = use_z; + + // Derived + GammaIdealFluid = eos_th.gamma; +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +c2p_1DEntropy::get_Ssq_Exact(const vec &mom, + const smat &gup) const { + + CCTK_REAL Ssq; + /* calculate S_squared */ + Ssq = mom(X) * (gup(X, X) * mom(X) + gup(X, Y) * mom(Y) + gup(X, Z) * mom(Z)); + Ssq += + mom(Y) * (gup(X, Y) * mom(X) + gup(Y, Y) * mom(Y) + gup(Y, Z) * mom(Z)); + Ssq += + mom(Z) * (gup(X, Z) * mom(X) + gup(Y, Z) * mom(Y) + gup(Z, Z) * mom(Z)); + + return Ssq; +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +c2p_1DEntropy::get_Bsq_Exact(const vec &B_up, + const smat &glo) const { + + vec B_low = calc_contraction(glo, B_up); + return calc_contraction(B_low, B_up); // Bsq +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +c2p_1DEntropy::get_BiSi_Exact(const vec &Bvec, + const vec &mom) const { + + return Bvec(X) * mom(X) + Bvec(Y) * mom(Y) + Bvec(Z) * mom(Z); // BiSi +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec +c2p_1DEntropy::get_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &v_up, + const smat &glo) const { + vec v_low = calc_contraction(glo, v_up); + CCTK_REAL vsq = calc_contraction(v_low, v_up); + CCTK_REAL VdotB = calc_contraction(v_low, B_up); + CCTK_REAL VdotBsq = VdotB * VdotB; + CCTK_REAL Bsq = get_Bsq_Exact(B_up, glo); + + CCTK_REAL w_lor = 1. / sqrt(1. - vsq); + CCTK_REAL bsq = ((Bsq) / (w_lor * w_lor)) + VdotBsq; + vec w_vsq_bsq{ w_lor, vsq, bsq }; + + return w_vsq_bsq; //{w_lor, vsq, bsq} +} + +/* Called by 1DEntropy */ + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +c2p_1DEntropy::xEntropyToPrim(CCTK_REAL xEntropy_Sol, CCTK_REAL Ssq, + CCTK_REAL Bsq, CCTK_REAL BiSi, const EOSType &eos_th, + prim_vars &pv, const cons_vars &cv, + const smat &gup, + const smat &glo) const { + // Density, entropy, Ye + pv.rho = xEntropy_Sol; + pv.entropy = cv.DEnt / cv.dens; + pv.Ye = cv.dYe / cv.dens; + + // Lorentz factor + pv.w_lor = cv.dens / xEntropy_Sol; + + // Pressure and epsilon + pv.press = + eos_th.press_from_valid_rho_kappa_ye(xEntropy_Sol, pv.entropy, pv.Ye); + pv.eps = eos_th.eps_from_valid_rho_kappa_ye(xEntropy_Sol, pv.entropy, pv.Ye); + + // Taken from WZ2Prim (2DNRNoble) + // Z_Sol = rho * h * w_lor * w_lor + CCTK_REAL Z_Sol = + (xEntropy_Sol + pv.eps * xEntropy_Sol + pv.press) * pv.w_lor * pv.w_lor; + + // TODO: Debug code to capture v>1, + // remove soon + if (use_zprim) { + + CCTK_REAL zx = + pv.w_lor * (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zx += pv.w_lor * BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + + CCTK_REAL zy = + pv.w_lor * (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zy += pv.w_lor * BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + + CCTK_REAL zz = + pv.w_lor * (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zz += pv.w_lor * BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + + CCTK_REAL zx_down = glo(X, X) * zx + glo(X, Y) * zy + glo(X, Z) * zz; + CCTK_REAL zy_down = glo(X, Y) * zx + glo(Y, Y) * zy + glo(Y, Z) * zz; + CCTK_REAL zz_down = glo(X, Z) * zx + glo(Y, Z) * zy + glo(Z, Z) * zz; + + CCTK_REAL Zsq = zx*zx_down + zy*zy_down + zz*zz_down; + + CCTK_REAL SafeLor = sqrt(1.0+Zsq); + + pv.vel(X) = zx/SafeLor; + pv.vel(Y) = zy/SafeLor; + pv.vel(Z) = zz/SafeLor; + + pv.w_lor = SafeLor; + + } else { + + pv.vel(X) = + (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(X) += BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Y) = + (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Y) += BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Z) = + (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Z) += BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + + } + + pv.Bvec = cv.dBvec; + + const vec Elow = calc_cross_product(pv.Bvec, pv.vel); + pv.E = calc_contraction(gup, Elow); +} + +// See Appendix A.4 of https://arxiv.org/pdf/1112.0568 +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +c2p_1DEntropy::funcRoot_1DEntropy(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, + CCTK_REAL x, const EOSType &eos_th, + const cons_vars &cv) const { + + // We already divided dens, DEnt and + // dYe by sqrt(gamma) + + const CCTK_REAL ent_loc = cv.DEnt / cv.dens; + const CCTK_REAL ye_loc = cv.dYe / cv.dens; + + // Compute h using entropy + const CCTK_REAL press_loc = + eos_th.press_from_valid_rho_kappa_ye(x, ent_loc, ye_loc); + + const CCTK_REAL eps_loc = + eos_th.eps_from_valid_rho_kappa_ye(x, ent_loc, ye_loc); + + // Compute (A60) using + // W = rho*h*lorentz*lorentz + const CCTK_REAL lor_loc = cv.dens / x; + const CCTK_REAL h_loc = (1.0 + eps_loc + press_loc / x); + const CCTK_REAL W = x * h_loc * lor_loc * lor_loc; + + // Compute (A61) + const CCTK_REAL Sfsq = + (W * W * Ssq + BiSi * BiSi * (Bsq + 2.0 * W)) / ((W + Bsq) * (W + Bsq)); + + // Compute (A62) + const CCTK_REAL rho = + cv.dens / sqrt(1.0 + Sfsq / (cv.dens * cv.dens * h_loc * h_loc)); + + return x - rho; +} + +template +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +c2p_1DEntropy::solve(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, + const smat &glo, c2p_report &rep) const { + + //ROOTSTAT status = ROOTSTAT::SUCCESS; + rep.iters = 0; + rep.adjust_cons = false; + rep.set_atmo = false; + rep.status = c2p_report::SUCCESS; + + /* Check validity of the 3-metric and compute its inverse */ + const CCTK_REAL spatial_detg = calc_det(glo); + const CCTK_REAL sqrt_detg = sqrt(spatial_detg); + + // if ((!isfinite(sqrt_detg)) || (sqrt_detg <= 0)) { + // rep.set_invalid_detg(sqrt_detg); + // set_to_nan(pv, cv); + // return; + //} + + // Check positive definiteness of spatial metric + // Sylvester's criterion, see + // https://en.wikipedia.org/wiki/Sylvester%27s_criterion + const bool minor1{ glo(X, X) > 0.0 }; + const bool minor2{ glo(X, X) * glo(Y, Y) - glo(X, Y) * glo(X, Y) > 0.0 }; + const bool minor3{ spatial_detg > 0.0 }; + + if (!(minor1 && minor2 && minor3)) { + rep.set_invalid_detg(sqrt_detg); + set_to_nan(pv, cv); + return; + } + + const smat gup = calc_inv(glo, spatial_detg); + + /* Undensitize the conserved vars */ + cv.dens /= sqrt_detg; + cv.tau /= sqrt_detg; + cv.mom /= sqrt_detg; + cv.dBvec /= sqrt_detg; + cv.dYe /= sqrt_detg; + cv.DEnt /= sqrt_detg; + + if (cv.dens <= atmo.rho_cut) { + rep.set_atmo_set(); + atmo.set(pv, cv, glo); + return; + } + + const CCTK_REAL Ssq = get_Ssq_Exact(cv.mom, gup); + const CCTK_REAL Bsq = get_Bsq_Exact(cv.dBvec, glo); + const CCTK_REAL BiSi = get_BiSi_Exact(cv.dBvec, cv.mom); + + // TODO: Is this check really necessary? + if ( (!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || + (!isfinite(BiSi)) || (!isfinite(cv.dYe)) || (!isfinite(cv.DEnt)) ) { + rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.dYe); + set_to_nan(pv, cv); + return; + } + + if (Bsq > Bsq_lim) { + rep.set_B_limit(Bsq); + set_to_nan(pv, cv); + return; + } + + // See Appendix A.4 of https://arxiv.org/pdf/1112.0568 + // Compute (A59) + CCTK_REAL a = cv.dens / sqrt(1.0 + Ssq / (cv.dens * cv.dens)); + CCTK_REAL b = cv.dens; + auto fn = [&](auto x) { + return funcRoot_1DEntropy(Ssq, Bsq, BiSi, x, eos_th, cv); + }; + + // Important! + // Algo::brent terminates if the following accuracy is achieved + // abs(x - y) <= eps * min(abs(x), abs(y)), + // where x and y are the values of the bracket and + // eps = std::ldexp(1, -minbits) = 1 * 2^{-minbits} + // This should probably be changed in Algo::brent + + // We want to set the tolerance to its correct parameter + const CCTK_REAL log2 = std::log(2.0); + const CCTK_INT minbits = int(abs(std::log(tolerance)) / log2); + const CCTK_REAL tolerance_0 = std::ldexp(double(1.0), -minbits); + // Old code: + // const CCTK_INT minbits = std::numeric_limits::digits - 4; + const CCTK_INT maxiters = maxIterations; + + auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); + + CCTK_REAL xEntropy_Sol = 0.5 * (result.first + result.second); + + xEntropyToPrim(xEntropy_Sol, Ssq, Bsq, BiSi, eos_th, pv, cv, gup, glo); + + // Check solution and calculate primitives + // if (rep.iters < maxiters && abs(fn(xEntropy_Sol)) < tolerance) { + /* + if (abs(result.first - result.second) <= + tolerance_0 * min(abs(result.first), abs(result.second))) { + rep.status = c2p_report::SUCCESS; + status = ROOTSTAT::SUCCESS; + } else { + // set status to root not converged + rep.set_root_conv(); + status = ROOTSTAT::NOT_CONVERGED; + } + */ + if (abs(result.first - result.second) > + tolerance_0 * min(abs(result.first), abs(result.second))) { + + // check primitives against conservatives + cons_vars cv_check; + cv_check.from_prim(pv, glo); + /* Undensitize the conserved vars */ + cv_check.dens /= sqrt_detg; + cv_check.tau /= sqrt_detg; + cv_check.mom /= sqrt_detg; + cv_check.dBvec /= sqrt_detg; + cv_check.dYe /= sqrt_detg; + cv_check.DEnt /= sqrt_detg; + + CCTK_REAL small = 1e-50; + + CCTK_REAL max_error = sqrt(max({pow((cv_check.dens-cv.dens)/(cv.dens+small),2.0), + pow((cv_check.mom(0)-cv.mom(0))/(cv.mom(0)+small),2.0), + pow((cv_check.mom(1)-cv.mom(1))/(cv.mom(1)+small),2.0), + pow((cv_check.mom(2)-cv.mom(2))/(cv.mom(2)+small),2.0), + pow((cv_check.tau-cv.tau)/(cv.tau+small),2.0)})); + + // reject only if mismatch in conservatives is inappropriate, else accept + if (max_error >= cons_error) { + // set status to root not converged + rep.set_root_conv(); + //status = ROOTSTAT::NOT_CONVERGED; + return; + } + } + + // Lower velocity + const vec v_low = calc_contraction(glo, pv.vel); + // Computing b^t : this is b^0 * alp + const CCTK_REAL bst = pv.w_lor * calc_contraction(pv.Bvec, v_low); + // Computing b^mu b_mu + const CCTK_REAL bs2 = (Bsq + bst * bst) / (pv.w_lor * pv.w_lor); + // Recompute tau + cv.tau = (pv.w_lor * pv.w_lor * (pv.rho * (1.0 + pv.eps) + pv.press + bs2) - + (pv.press + 0.5 * bs2) - bst * bst) - + cv.dens; + + // set to atmo if computed rho is below floor density + if (pv.rho < atmo.rho_cut) { + rep.set_atmo_set(); + atmo.set(pv, cv, glo); + return; + } + + c2p::prims_floors_and_ceilings(eos_th,pv,cv,glo,rep); + + // Recompute cons if prims have been adjusted + if (rep.adjust_cons) { + cv.from_prim(pv, glo); + } else { + /* Densitize the conserved vars again*/ + cv.dens *= sqrt_detg; + cv.tau *= sqrt_detg; + cv.mom *= sqrt_detg; + cv.dBvec *= sqrt_detg; + cv.dYe *= sqrt_detg; + cv.DEnt *= sqrt_detg; + } + +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void +c2p_1DEntropy::set_to_nan(prim_vars &pv, cons_vars &cv) const { + pv.set_to_nan(); + cv.set_to_nan(); +} + +/* Destructor */ +CCTK_HOST CCTK_DEVICE +CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DEntropy::~c2p_1DEntropy() { + // How to destruct properly a vector? +} +} // namespace Con2PrimFactory + +#endif diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 5e1f7c43..db53692c 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -14,20 +14,22 @@ public: /* Constructor */ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, - CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len); + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z); CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Ssq_Exact(vec &mom, + get_Ssq_Exact(const vec &mom, const smat &gup) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Bsq_Exact(vec &B_up, + get_Bsq_Exact(const vec &B_up, const smat &glo) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_BiSi_Exact(vec &Bvec, vec &mom) const; + get_BiSi_Exact(const vec &Bvec, const vec &mom) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec - get_WLorentz_vsq_bsq_Seeds(vec &B_up, - vec &v_up, + get_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &v_up, const smat &glo) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void set_to_nan(prim_vars &pv, cons_vars &cv) const; @@ -35,18 +37,18 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq, - CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv, - cons_vars cv, const smat &gup, + CCTK_REAL BiSi, const EOSType &eos_th, prim_vars &pv, + const cons_vars &cv, const smat &gup, const smat &glo) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, - CCTK_REAL x, EOSType &eos_th, cons_vars &cv) const; + CCTK_REAL x, const EOSType &eos_th, const cons_vars &cv) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars &cv, + solve(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, const smat &glo, c2p_report &rep) const; /* Destructor */ @@ -57,23 +59,33 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela::c2p_1DPalenzuela( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, - CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) { + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z) { - GammaIdealFluid = eos_th.gamma; + // Base + atmo = atm; maxIterations = maxIter; tolerance = tol; - rho_strict = rho_str; - ye_lenient = ye_len; + alp_thresh = alp_thresh_in; + cons_error = consError; vw_lim = vwlim; w_lim = sqrt(1.0 + vw_lim * vw_lim); v_lim = vw_lim / w_lim; Bsq_lim = B_lim * B_lim; - atmo = atm; + rho_BH = rho_BH_in; + eps_BH = eps_BH_in; + vwlim_BH = vwlim_BH_in; + ye_lenient = ye_len; + use_zprim = use_z; + + // Derived + GammaIdealFluid = eos_th.gamma; } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_1DPalenzuela::get_Ssq_Exact(vec &mom, +c2p_1DPalenzuela::get_Ssq_Exact(const vec &mom, const smat &gup) const { CCTK_REAL Ssq; @@ -84,14 +96,11 @@ c2p_1DPalenzuela::get_Ssq_Exact(vec &mom, Ssq += mom(Z) * (gup(X, Z) * mom(X) + gup(Y, Z) * mom(Y) + gup(Z, Z) * mom(Z)); - if ((Ssq < 0.) && (fabs(Ssq) < 1.0e-13)) { - Ssq = fabs(Ssq); - } return Ssq; } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_1DPalenzuela::get_Bsq_Exact(vec &B_up, +c2p_1DPalenzuela::get_Bsq_Exact(const vec &B_up, const smat &glo) const { vec B_low = calc_contraction(glo, B_up); @@ -99,15 +108,15 @@ c2p_1DPalenzuela::get_Bsq_Exact(vec &B_up, } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_1DPalenzuela::get_BiSi_Exact(vec &Bvec, - vec &mom) const { +c2p_1DPalenzuela::get_BiSi_Exact(const vec &Bvec, + const vec &mom) const { return Bvec(X) * mom(X) + Bvec(Y) * mom(Y) + Bvec(Z) * mom(Z); // BiSi } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec c2p_1DPalenzuela::get_WLorentz_vsq_bsq_Seeds( - vec &B_up, vec &v_up, + const vec &B_up, const vec &v_up, const smat &glo) const { vec v_low = calc_contraction(glo, v_up); CCTK_REAL vsq = calc_contraction(v_low, v_up); @@ -115,10 +124,6 @@ c2p_1DPalenzuela::get_WLorentz_vsq_bsq_Seeds( CCTK_REAL VdotBsq = VdotB * VdotB; CCTK_REAL Bsq = get_Bsq_Exact(B_up, glo); - if ((vsq < 0.) && (fabs(vsq) < 1.0e-13)) { - vsq = fabs(vsq); - } - CCTK_REAL w_lor = 1. / sqrt(1. - vsq); CCTK_REAL bsq = ((Bsq) / (w_lor * w_lor)) + VdotBsq; vec w_vsq_bsq{w_lor, vsq, bsq}; @@ -132,8 +137,8 @@ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, - EOSType &eos_th, prim_vars &pv, - cons_vars cv, const smat &gup, + const EOSType &eos_th, prim_vars &pv, + const cons_vars &cv, const smat &gup, const smat &glo) const { const CCTK_REAL qPalenzuela = cv.tau / cv.dens; const CCTK_REAL rPalenzuela = Ssq / pow(cv.dens, 2); @@ -166,27 +171,66 @@ c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, // Taken from WZ2Prim (2DNRNoble) CCTK_REAL Z_Sol = xPalenzuela_Sol * pv.rho * W_sol; - pv.vel(X) = - (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(X) += BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + // TODO: Debug code to capture v>1, + // remove soon + if (use_zprim) { + + CCTK_REAL zx = + W_sol * (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zx += W_sol * BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + + CCTK_REAL zy = + W_sol * (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zy += W_sol * BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + + CCTK_REAL zz = + W_sol * (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zz += W_sol * BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); - pv.vel(Y) = - (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(Y) += BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + CCTK_REAL zx_down = glo(X, X) * zx + glo(X, Y) * zy + glo(X, Z) * zz; + CCTK_REAL zy_down = glo(X, Y) * zx + glo(Y, Y) * zy + glo(Y, Z) * zz; + CCTK_REAL zz_down = glo(X, Z) * zx + glo(Y, Z) * zy + glo(Z, Z) * zz; - pv.vel(Z) = - (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(Z) += BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + CCTK_REAL Zsq = zx*zx_down + zy*zy_down + zz*zz_down; - pv.w_lor = W_sol; + CCTK_REAL SafeLor = sqrt(1.0+Zsq); + + pv.vel(X) = zx/SafeLor; + pv.vel(Y) = zy/SafeLor; + pv.vel(Z) = zz/SafeLor; + + pv.w_lor = SafeLor; + + } else { + + pv.vel(X) = + (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(X) += BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Y) = + (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Y) += BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Z) = + (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Z) += BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + + pv.w_lor = W_sol; + + } pv.Ye = cv.dYe / cv.dens; pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_th.kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.Bvec = cv.dBvec; const vec Elow = calc_cross_product(pv.Bvec, pv.vel); @@ -197,7 +241,7 @@ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi, CCTK_REAL x, - EOSType &eos_th, cons_vars &cv) const { + const EOSType &eos_th, const cons_vars &cv) const { // computes f(x) from x and q,r,s,t const CCTK_REAL qPalenzuela = cv.tau / cv.dens; const CCTK_REAL rPalenzuela = Ssq / pow(cv.dens, 2); @@ -231,11 +275,11 @@ c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, +c2p_1DPalenzuela::solve(const EOSType &eos_th, prim_vars &pv, cons_vars &cv, const smat &glo, c2p_report &rep) const { - ROOTSTAT status = ROOTSTAT::SUCCESS; + //ROOTSTAT status = ROOTSTAT::SUCCESS; rep.iters = 0; rep.adjust_cons = false; rep.set_atmo = false; @@ -244,11 +288,26 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, /* Check validity of the 3-metric and compute its inverse */ const CCTK_REAL spatial_detg = calc_det(glo); const CCTK_REAL sqrt_detg = sqrt(spatial_detg); - if ((!isfinite(sqrt_detg)) || (sqrt_detg <= 0)) { + + //if ((!isfinite(sqrt_detg)) || (sqrt_detg <= 0)) { + // rep.set_invalid_detg(sqrt_detg); + // set_to_nan(pv, cv); + // return; + //} + + // Check positive definiteness of spatial metric + // Sylvester's criterion, see + // https://en.wikipedia.org/wiki/Sylvester%27s_criterion + const bool minor1{ glo(X, X) > 0.0 }; + const bool minor2{ glo(X, X) * glo(Y, Y) - glo(X, Y) * glo(X, Y) > 0.0 }; + const bool minor3{ spatial_detg > 0.0 }; + + if (!(minor1 && minor2 && minor3)) { rep.set_invalid_detg(sqrt_detg); set_to_nan(pv, cv); return; } + const smat gup = calc_inv(glo, spatial_detg); /* Undensitize the conserved vars */ @@ -257,6 +316,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cv.mom /= sqrt_detg; cv.dBvec /= sqrt_detg; cv.dYe /= sqrt_detg; + cv.DEnt /= sqrt_detg; if (cv.dens <= atmo.rho_cut) { rep.set_atmo_set(); @@ -264,28 +324,14 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, return; } - // compute primitive B seed from conserved B of current time step for better - // guess - pv_seeds.Bvec = cv.dBvec; - const CCTK_REAL Ssq = get_Ssq_Exact(cv.mom, gup); - const CCTK_REAL Bsq = get_Bsq_Exact(pv_seeds.Bvec, glo); - const CCTK_REAL BiSi = get_BiSi_Exact(pv_seeds.Bvec, cv.mom); - const vec w_vsq_bsq = get_WLorentz_vsq_bsq_Seeds( - pv_seeds.Bvec, pv_seeds.vel, glo); // this also recomputes pv_seeds.w_lor - pv_seeds.w_lor = w_vsq_bsq(0); - // CCTK_REAL vsq_seed = w_vsq_bsq(1); - // const CCTK_REAL bsq = w_vsq_bsq(2); - - if ((!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || - (!isfinite(BiSi)) || (!isfinite(cv.dYe))) { - rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.dYe); - set_to_nan(pv, cv); - return; - } + const CCTK_REAL Bsq = get_Bsq_Exact(cv.dBvec, glo); + const CCTK_REAL BiSi = get_BiSi_Exact(cv.dBvec,cv.mom); - if (Bsq < 0) { - rep.set_neg_Bsq(Bsq); + // TODO: Is this check really necessary? + if ( (!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || + (!isfinite(BiSi)) || (!isfinite(cv.dYe)) || (!isfinite(cv.DEnt)) ) { + rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.dYe); set_to_nan(pv, cv); return; } @@ -296,13 +342,23 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, return; } - /* update rho seed from cv and wlor */ - // rho consistent with cv.rho should be better guess than rho from last - // timestep - pv_seeds.rho = cv.dens / pv_seeds.w_lor; - // Find x, this is the recovery process - const CCTK_INT minbits = std::numeric_limits::digits - 4; + //const CCTK_INT minbits = std::numeric_limits::digits - 4; + //const CCTK_INT maxiters = maxIterations; + + // Important! + // Algo::brent terminates if the following accuracy is achieved + // abs(x - y) <= eps * min(abs(x), abs(y)), + // where x and y are the values of the bracket and + // eps = std::ldexp(1, -minbits) = 1 * 2^{-minbits} + // This should probably be changed in Algo::brent + + // We want to set the tolerance to its correct parameter + const CCTK_REAL log2 = std::log(2.0); + const CCTK_INT minbits = int(abs(std::log(tolerance)) / log2); + const CCTK_REAL tolerance_0 = std::ldexp(double(1.0), -minbits); + // Old code: + // const CCTK_INT minbits = std::numeric_limits::digits - 4; const CCTK_INT maxiters = maxIterations; CCTK_REAL qPalenzuela = cv.tau / cv.dens; @@ -316,27 +372,74 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, }; auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters); - // Pick best solution - CCTK_REAL xPalenzuela_Sol; - if (abs(fn(result.first)) < abs(fn(result.second))) { - xPalenzuela_Sol = result.first; - } else { - xPalenzuela_Sol = result.second; - } + CCTK_REAL xPalenzuela_Sol = 0.5 * (result.first + result.second); + + //if (abs(fn(result.first)) < abs(fn(result.second))) { + // xPalenzuela_Sol = result.first; + //} else { + // xPalenzuela_Sol = result.second; + //} // Check solution and calculate primitives // TODO:check if to pass result.first or xPalenzuela_Sol - if (rep.iters < maxiters && abs(fn(xPalenzuela_Sol)) < tolerance) { - rep.status = c2p_report::SUCCESS; - status = ROOTSTAT::SUCCESS; - } else { + //if (rep.iters < maxiters && abs(fn(xPalenzuela_Sol)) < tolerance) { + // rep.status = c2p_report::SUCCESS; + // status = ROOTSTAT::SUCCESS; + //} else { // set status to root not converged - rep.set_root_conv(); - status = ROOTSTAT::NOT_CONVERGED; - } + // rep.set_root_conv(); + // status = ROOTSTAT::NOT_CONVERGED; + //} xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_th, pv, cv, gup, glo); + + // General comment: + // One could think of expressing the following condition + // in a way that is "safe" against NaNs and infs. First, this only makes + // sense if we want these values to be considered as failures which should + // be treated as "not converged". + // + // inf: Since inf behaves like a large valid number nothing special needs + // to be done except of rewriting the argument of the if condition such that + // possible infs are present only on one side of the comparison, eg + // abs(difference)/abs(normalization) > tolerance_0 + // + // NaN: If the argument of if (...) is NaN, it usually evaluates to false. + // Here, we would need to rewrite the logic a little bit. + if (abs(result.first - result.second) > + tolerance_0 * min(abs(result.first), abs(result.second))) { + + // check primitives against conservatives + cons_vars cv_check; + cv_check.from_prim(pv, glo); + + /* Undensitize the conserved vars */ + cv_check.dens /= sqrt_detg; + cv_check.tau /= sqrt_detg; + cv_check.mom /= sqrt_detg; + cv_check.dBvec /= sqrt_detg; + cv_check.dYe /= sqrt_detg; + cv_check.DEnt /= sqrt_detg; + + CCTK_REAL small = 1e-50; + + CCTK_REAL max_error = sqrt(max({pow((cv_check.dens-cv.dens)/(cv.dens+small),2.0), + pow((cv_check.mom(0)-cv.mom(0))/(cv.mom(0)+small),2.0), + pow((cv_check.mom(1)-cv.mom(1))/(cv.mom(1)+small),2.0), + pow((cv_check.mom(2)-cv.mom(2))/(cv.mom(2)+small),2.0), + pow((cv_check.tau-cv.tau)/(cv.tau+small),2.0)})); + + // reject only if mismatch in conservatives is inappropriate, else accept + if (max_error >= cons_error) { + + // set status to root not converged + rep.set_root_conv(); + //status = ROOTSTAT::NOT_CONVERGED; + return; + } + } + // set to atmo if computed rho is below floor density if (pv.rho < atmo.rho_cut) { rep.set_atmo_set(); @@ -344,66 +447,19 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, return; } - // check the validity of the computed eps - auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); - if (pv.eps > rgeps.max) { - //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n"); - rep.adjust_cons = true; - if (pv.rho >= rho_strict) { - rep.set_range_eps(pv.eps); // sets adjust_cons to false by default - rep.adjust_cons = true; - set_to_nan(pv, cv); - return; - } - } else if (pv.eps < rgeps.min) { - - /* - printf( - "(pv.eps < rgeps.min) is true! pv.eps, rgeps.min: %26.16e, %26.16e \n", - pv.eps, rgeps.min); - printf("Adjusting cons.. \n"); - */ - pv.eps = rgeps.min; - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - rep.set_range_eps(rgeps.min); // sets adjust_cons to true - } - - // TODO: check validity for Ye - - // check if computed velocities are within the specified limit - vec v_low = calc_contraction(glo, pv.vel); - CCTK_REAL vsq_Sol = calc_contraction(v_low, pv.vel); - CCTK_REAL sol_v = sqrt(vsq_Sol); - if (sol_v > v_lim) { - /* - printf("(sol_v > v_lim) is true! \n"); - printf("sol_v, v_lim: %26.16e, %26.16e \n", sol_v, v_lim); - */ - pv.rho = cv.dens / w_lim; - if (pv.rho >= rho_strict) { - rep.set_speed_limit({sol_v, sol_v, sol_v}); - set_to_nan(pv, cv); - return; - } - pv.vel *= v_lim / sol_v; - pv.w_lor = w_lim; - pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max); - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - - rep.adjust_cons = true; - } + c2p::prims_floors_and_ceilings(eos_th,pv,cv,glo,rep); // Recompute cons if prims have been adjusted if (rep.adjust_cons) { cv.from_prim(pv, glo); - } - else { - /* If not adjusted, densitize back the original conserved vars */ + } else { + /* Densitize the conserved vars again*/ cv.dens *= sqrt_detg; cv.tau *= sqrt_detg; cv.mom *= sqrt_detg; cv.dBvec *= sqrt_detg; cv.dYe *= sqrt_detg; + cv.DEnt *= sqrt_detg; } } diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx index 0185392b..b3f2a694 100644 --- a/Con2PrimFactory/src/c2p_2DNoble.hxx +++ b/Con2PrimFactory/src/c2p_2DNoble.hxx @@ -11,24 +11,34 @@ class c2p_2DNoble : public c2p { public: /* Some attributes */ CCTK_REAL GammaIdealFluid; + CCTK_REAL Zmin; /* Constructor */ template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, - CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len); + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z); + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Ssq_Exact(vec &mom, + get_Ssq_Exact(const vec &mom, const smat &gup) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Bsq_Exact(vec &B_up, + get_Bsq_Exact(const vec &B_up, const smat &glo) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_BiSi_Exact(vec &Bvec, vec &mom) const; + get_BiSi_Exact(const vec &Bvec, const vec &mom) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec - get_WLorentz_vsq_bsq_Seeds(vec &B_up, - vec &v_up, + get_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &v_up, const smat &glo) const; + // TODO: Debug function to capture v>1, + // remove soon + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec + getZ_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &z_up, + const smat &glo) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void set_to_nan(prim_vars &pv, cons_vars &cv) const; @@ -37,19 +47,20 @@ public: get_Z_Seed(CCTK_REAL rho, CCTK_REAL eps, CCTK_REAL press, CCTK_REAL w_lor) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_Press_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, cons_vars &cv) const; + get_Press_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, const cons_vars &cv) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL get_dPdZ_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq) const; CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL - get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, cons_vars &cv) const; + get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, const cons_vars &cv) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, CCTK_REAL BiSi, - EOSType &eos_th, prim_vars &pv, cons_vars cv, - const smat &gup) const; + const EOSType &eos_th, prim_vars &pv, const cons_vars &cv, + const smat &gup, + const smat &glo) const; template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars &cv, + solve(const EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars &cv, const smat &glo, c2p_report &rep) const; /* Destructor */ @@ -60,49 +71,57 @@ public: template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble::c2p_2DNoble( - EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, - CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) { + const EOSType &eos_th, const atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol, + CCTK_REAL alp_thresh_in, CCTK_REAL consError, + CCTK_REAL vwlim, CCTK_REAL B_lim, CCTK_REAL rho_BH_in, CCTK_REAL eps_BH_in, CCTK_REAL vwlim_BH_in, + bool ye_len, bool use_z) { - GammaIdealFluid = eos_th.gamma; + // Base + atmo = atm; maxIterations = maxIter; tolerance = tol; - rho_strict = rho_str; - ye_lenient = ye_len; + alp_thresh = alp_thresh_in; + cons_error = consError; vw_lim = vwlim; w_lim = sqrt(1.0 + vw_lim * vw_lim); v_lim = vw_lim / w_lim; Bsq_lim = B_lim * B_lim; - atmo = atm; + rho_BH = rho_BH_in; + eps_BH = eps_BH_in; + vwlim_BH = vwlim_BH_in; + ye_lenient = ye_len; + use_zprim = use_z; + + // Derived + GammaIdealFluid = eos_th.gamma; + Zmin = eos_th.rgrho.min; } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_2DNoble::get_Ssq_Exact(vec &mom_low, +c2p_2DNoble::get_Ssq_Exact(const vec &mom_low, const smat &gup) const { vec mom_up = calc_contraction(gup, mom_low); CCTK_REAL Ssq = calc_contraction(mom_low, mom_up); - if ((Ssq < 0.) && (fabs(Ssq) < 1.0e-13)) { - Ssq = fabs(Ssq); - } return Ssq; } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_2DNoble::get_Bsq_Exact(vec &B_up, +c2p_2DNoble::get_Bsq_Exact(const vec &B_up, const smat &glo) const { vec B_low = calc_contraction(glo, B_up); return calc_contraction(B_low, B_up); // Bsq } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL -c2p_2DNoble::get_BiSi_Exact(vec &Bvec, - vec &mom) const { +c2p_2DNoble::get_BiSi_Exact(const vec &Bvec, + const vec &mom) const { return calc_contraction(mom, Bvec); // BiS^i } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -c2p_2DNoble::get_WLorentz_vsq_bsq_Seeds(vec &B_up, - vec &v_up, +c2p_2DNoble::get_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &v_up, const smat &glo) const { vec v_low = calc_contraction(glo, v_up); CCTK_REAL vsq = calc_contraction(v_low, v_up); @@ -110,10 +129,6 @@ c2p_2DNoble::get_WLorentz_vsq_bsq_Seeds(vec &B_up, CCTK_REAL VdotBsq = VdotB * VdotB; CCTK_REAL Bsq = get_Bsq_Exact(B_up, glo); - if ((vsq < 0.) && (fabs(vsq) < 1.0e-13)) { - vsq = fabs(vsq); - } - CCTK_REAL w_lor = 1. / sqrt(1. - vsq); CCTK_REAL bsq = ((Bsq) / (w_lor * w_lor)) + VdotBsq; vec w_vsq_bsq{w_lor, vsq, bsq}; @@ -121,6 +136,28 @@ c2p_2DNoble::get_WLorentz_vsq_bsq_Seeds(vec &B_up, return w_vsq_bsq; //{w_lor, vsq, bsq} } +// TODO: Debug function to capture v>1, +// remove soon +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec +c2p_2DNoble::getZ_WLorentz_vsq_bsq_Seeds(const vec &B_up, + const vec &z_up, + const smat &glo) const { + vec z_low = calc_contraction(glo, z_up); + CCTK_REAL zsq = calc_contraction(z_low, z_up); + + CCTK_REAL w_lor = sqrt(1. + zsq); + CCTK_REAL vsq = min(zsq/w_lor/w_lor,1.-1.e-15); + + CCTK_REAL VdotB = calc_contraction(z_low, B_up)/w_lor; + CCTK_REAL VdotBsq = VdotB * VdotB; + CCTK_REAL Bsq = get_Bsq_Exact(B_up, glo); + + CCTK_REAL bsq = ((Bsq) / (w_lor * w_lor)) + VdotBsq; + vec w_vsq_bsq{w_lor, vsq, bsq}; + + return w_vsq_bsq; //{w_lor, vsq, bsq} +} + /* Called by 2dNRNoble */ CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL @@ -131,7 +168,7 @@ c2p_2DNoble::get_Z_Seed(CCTK_REAL rho, CCTK_REAL eps, CCTK_REAL press, CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL c2p_2DNoble::get_Press_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, - cons_vars &cv) const { + const cons_vars &cv) const { return ((Z * (1.0 - Vsq) - cv.dens * sqrt(1.0 - Vsq)) * (GammaIdealFluid - 1.0) / (GammaIdealFluid)); } @@ -143,7 +180,7 @@ c2p_2DNoble::get_dPdZ_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq) const { CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL c2p_2DNoble::get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, - cons_vars &cv) const { + const cons_vars &cv) const { return ((-Z + cv.dens / (2.0 * sqrt(1.0 - Vsq))) * (GammaIdealFluid - 1.0) / GammaIdealFluid); } @@ -151,35 +188,76 @@ c2p_2DNoble::get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq, template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, - CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv, - cons_vars cv, const smat &gup) const { + CCTK_REAL BiSi, const EOSType &eos_th, prim_vars &pv, + const cons_vars &cv, const smat &gup, + const smat &glo) const { CCTK_REAL W_Sol = 1.0 / sqrt(1.0 - vsq_Sol); pv.rho = cv.dens / W_Sol; - pv.vel(X) = - (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(X) += BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + // TODO: Debug code to capture v>1, + // remove soon + if (use_zprim) { + + CCTK_REAL zx = + W_Sol * (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zx += W_Sol * BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); - pv.vel(Y) = - (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(Y) += BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + CCTK_REAL zy = + W_Sol * (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zy += W_Sol * BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); - pv.vel(Z) = - (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / - (Z_Sol + Bsq); - pv.vel(Z) += BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + CCTK_REAL zz = + W_Sol * (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + zz += W_Sol * BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); - pv.w_lor = W_Sol; + CCTK_REAL zx_down = glo(X, X) * zx + glo(X, Y) * zy + glo(X, Z) * zz; + CCTK_REAL zy_down = glo(X, Y) * zx + glo(Y, Y) * zy + glo(Y, Z) * zz; + CCTK_REAL zz_down = glo(X, Z) * zx + glo(Y, Z) * zy + glo(Z, Z) * zz; - pv.eps = (Z_Sol * (1. - vsq_Sol) / pv.rho - 1.0) / GammaIdealFluid; + CCTK_REAL Zsq = zx*zx_down + zy*zy_down + zz*zz_down; + + CCTK_REAL SafeLor = sqrt(1.0+Zsq); + + pv.vel(X) = zx/SafeLor; + pv.vel(Y) = zy/SafeLor; + pv.vel(Z) = zz/SafeLor; + + pv.w_lor = SafeLor; + + } else { + + pv.vel(X) = + (gup(X, X) * cv.mom(X) + gup(X, Y) * cv.mom(Y) + gup(X, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(X) += BiSi * cv.dBvec(X) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Y) = + (gup(X, Y) * cv.mom(X) + gup(Y, Y) * cv.mom(Y) + gup(Y, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Y) += BiSi * cv.dBvec(Y) / (Z_Sol * (Z_Sol + Bsq)); + + pv.vel(Z) = + (gup(X, Z) * cv.mom(X) + gup(Y, Z) * cv.mom(Y) + gup(Z, Z) * cv.mom(Z)) / + (Z_Sol + Bsq); + pv.vel(Z) += BiSi * cv.dBvec(Z) / (Z_Sol * (Z_Sol + Bsq)); + + pv.w_lor = W_Sol; + + } + + //pv.eps = (Z_Sol * (1. - vsq_Sol) / pv.rho - 1.0) / GammaIdealFluid; + pv.eps = (Z_Sol / pv.w_lor / pv.w_lor / pv.rho - 1.0) / GammaIdealFluid; pv.Ye = cv.dYe / cv.dens; pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.entropy = eos_th.kappa_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); + pv.Bvec = cv.dBvec; const vec Elow = calc_cross_product(pv.Bvec, pv.vel); @@ -188,11 +266,11 @@ c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, template CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void -c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, +c2p_2DNoble::solve(const EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars &cv, const smat &glo, c2p_report &rep) const { - ROOTSTAT status = ROOTSTAT::SUCCESS; + //ROOTSTAT status = ROOTSTAT::SUCCESS; rep.iters = 0; rep.adjust_cons = false; rep.set_atmo = false; @@ -202,11 +280,25 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, const CCTK_REAL spatial_detg = calc_det(glo); const CCTK_REAL sqrt_detg = sqrt(spatial_detg); - if ((!isfinite(sqrt_detg)) || (sqrt_detg <= 0)) { + //if ((!isfinite(sqrt_detg)) || (sqrt_detg <= 0)) { + // rep.set_invalid_detg(sqrt_detg); + // set_to_nan(pv, cv); + // return; + //} + + // Check positive definiteness of spatial metric + // Sylvester's criterion, see + // https://en.wikipedia.org/wiki/Sylvester%27s_criterion + const bool minor1{ glo(X, X) > 0.0 }; + const bool minor2{ glo(X, X) * glo(Y, Y) - glo(X, Y) * glo(X, Y) > 0.0 }; + const bool minor3{ spatial_detg > 0.0 }; + + if (!(minor1 && minor2 && minor3)) { rep.set_invalid_detg(sqrt_detg); set_to_nan(pv, cv); return; } + const smat gup = calc_inv(glo, spatial_detg); /* Undensitize the conserved vars */ @@ -215,6 +307,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cv.mom /= sqrt_detg; cv.dBvec /= sqrt_detg; cv.dYe /= sqrt_detg; + cv.DEnt /= sqrt_detg; if (cv.dens <= atmo.rho_cut) { rep.set_atmo_set(); @@ -229,21 +322,35 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, const CCTK_REAL Ssq = get_Ssq_Exact(cv.mom, gup); const CCTK_REAL Bsq = get_Bsq_Exact(pv_seeds.Bvec, glo); const CCTK_REAL BiSi = get_BiSi_Exact(pv_seeds.Bvec, cv.mom); - const vec w_vsq_bsq = get_WLorentz_vsq_bsq_Seeds( + + vec w_vsq_bsq; + CCTK_REAL vsq_seed; + + // TODO: Debug code to capture v>1, + // remove soon + if (use_zprim) { + + vec zvec = pv_seeds.vel*pv_seeds.w_lor; + + w_vsq_bsq = getZ_WLorentz_vsq_bsq_Seeds( + pv_seeds.Bvec, zvec, glo); // this also recomputes pv_seeds.w_lor + + pv_seeds.w_lor = w_vsq_bsq(0); + vsq_seed = w_vsq_bsq(1); + + } else { + + w_vsq_bsq = get_WLorentz_vsq_bsq_Seeds( pv_seeds.Bvec, pv_seeds.vel, glo); // this also recomputes pv_seeds.w_lor - pv_seeds.w_lor = w_vsq_bsq(0); - CCTK_REAL vsq_seed = w_vsq_bsq(1); - // const CCTK_REAL bsq = w_vsq_bsq(2); // small b^2 - if ((!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || - (!isfinite(BiSi)) || (!isfinite(cv.dYe))) { - rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.dYe); - set_to_nan(pv, cv); - return; + pv_seeds.w_lor = w_vsq_bsq(0); + vsq_seed = w_vsq_bsq(1); } - if (Bsq < 0) { - rep.set_neg_Bsq(Bsq); + // TODO: Is this check really necessary? + if ( (!isfinite(cv.dens)) || (!isfinite(Ssq)) || (!isfinite(Bsq)) || + (!isfinite(BiSi)) || (!isfinite(cv.dYe)) || (!isfinite(cv.DEnt)) ) { + rep.set_nans_in_cons(cv.dens, Ssq, Bsq, BiSi, cv.dYe); set_to_nan(pv, cv); return; } @@ -277,15 +384,15 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, x[0] = fabs(Z_Seed); x[1] = vsq_seed; - // printf("x[0], x[1]: %f, %f \n", x[0], x[1]); - /* initialize old values */ x_old[0] = x[0]; x_old[1] = x[1]; /* Start Recovery with 2D NR Solver */ - const CCTK_INT n = 2; - const CCTK_REAL dv = (1. - 1.e-15); + constexpr CCTK_INT n = 2; + constexpr CCTK_REAL dv = (1. - 1.e-10); + constexpr CCTK_REAL dw = 1. / (1. - dv); + CCTK_REAL dx[n]; CCTK_REAL fjac[n][n]; CCTK_REAL resid[n]; @@ -294,6 +401,26 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, CCTK_REAL df = 1.; CCTK_REAL f = 1.; + /* make sure that x[] is physical */ + if (x[1] < 0.0) { + x[1] = 0.0; + } + + else { + if (x[1] >= dv) { + x[0] *= dw*(1.0-x[1]); + x[1] = dv; + } + } + + if (x[0] < Zmin) { + x[0] = Zmin; + } else { + if (x[0] > 1e20) { + x[0] = x_old[0]; + } + } + CCTK_INT k; for (k = 1; k <= maxIterations; k++) { @@ -304,6 +431,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, const CCTK_REAL Z = x[0]; const CCTK_REAL invZ = 1.0 / Z; const CCTK_REAL Vsq = x[1]; + const CCTK_REAL Sdotn = -(cv.tau + cv.dens); const CCTK_REAL p_tmp = get_Press_funcZVsq(Z, Vsq, cv); const CCTK_REAL dPdvsq = get_dPdVsq_funcZVsq(Z, Vsq, cv); @@ -344,35 +472,48 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, errx = (x[0] == 0.) ? fabs(dx[0]) : fabs(dx[0] / x[0]); /* make sure that the new x[] is physical */ - if (x[0] < 0.0) { - x[0] = fabs(x[0]); - } else { - if (x[0] > 1e20) { - x[0] = x_old[0]; - } - } - if (x[1] < 0.0) { x[1] = 0.0; } else { - if (x[1] >= 1.0) { + if (x[1] >= dv) { + x[0] *= dw*(1.0-x[1]); x[1] = dv; } } + + if (x[0] < Zmin) { + x[0] = Zmin; + } else { + if (x[0] > 1e20) { + x[0] = x_old[0]; + } + } + + if (fabs(errx) <= tolerance) { + break; + } + } // storing number of iterations taken to find the root rep.iters = k; - if (fabs(errx) <= tolerance) { - rep.status = c2p_report::SUCCESS; - status = ROOTSTAT::SUCCESS; - } else { + //if (fabs(errx) <= tolerance) { + //rep.status = c2p_report::SUCCESS; + //status = ROOTSTAT::SUCCESS; + //} else { + // set status to root not converged + //rep.set_root_conv(); + //status = ROOTSTAT::NOT_CONVERGED; + //} + + if (fabs(errx) > tolerance) { // set status to root not converged rep.set_root_conv(); - status = ROOTSTAT::NOT_CONVERGED; + //status = ROOTSTAT::NOT_CONVERGED; + return; } // Check for bad untrapped divergences @@ -387,7 +528,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, CCTK_REAL vsq_Sol = x[1]; /* Write prims if C2P succeeded */ - WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_th, pv, cv, gup); + WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_th, pv, cv, gup, glo); // set to atmo if computed rho is below floor density if (pv.rho < atmo.rho_cut) { @@ -396,51 +537,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, return; } - // check the validity of the computed eps - auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye); - if (pv.eps > rgeps.max) { - //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n"); - rep.adjust_cons = true; - if (pv.rho >= rho_strict) { - rep.set_range_eps(pv.eps); // sets adjust_cons to false by default - rep.adjust_cons = true; - set_to_nan(pv, cv); - return; - } - } else if (pv.eps < rgeps.min) { - /* - printf( - "(pv.eps < rgeps.min) is true! pv.eps, rgeps.min: %26.16e, %26.16e \n", - pv.eps, rgeps.min); - printf(" Not adjusting cons.. \n"); - */ - pv.eps = rgeps.min; - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - rep.set_range_eps(rgeps.min); // sets adjust_cons to true by default - } - - // TODO: check validity for Ye - - // check if computed velocities are within the specified limit - CCTK_REAL sol_v = sqrt(vsq_Sol); - if (sol_v > v_lim) { - /* - printf("(sol_v > v_lim) is true! \n"); - printf("sol_v, v_lim: %26.16e, %26.16e \n", sol_v, v_lim); - */ - pv.rho = cv.dens / w_lim; - if (pv.rho >= rho_strict) { - rep.set_speed_limit({sol_v, sol_v, sol_v}); - set_to_nan(pv, cv); - return; - } - pv.vel *= v_lim / sol_v; - pv.w_lor = w_lim; - pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max); - pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye); - - rep.adjust_cons = true; - } + c2p::prims_floors_and_ceilings(eos_th,pv,cv,glo,rep); // Recompute cons if prims have been adjusted if (rep.adjust_cons) { @@ -453,6 +550,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cv.mom *= sqrt_detg; cv.dBvec *= sqrt_detg; cv.dYe *= sqrt_detg; + cv.DEnt *= sqrt_detg; } } diff --git a/Con2PrimFactory/src/cons.hxx b/Con2PrimFactory/src/cons.hxx index 9f0019bf..1067af04 100644 --- a/Con2PrimFactory/src/cons.hxx +++ b/Con2PrimFactory/src/cons.hxx @@ -17,6 +17,7 @@ struct cons_vars { vec mom; CCTK_REAL tau; CCTK_REAL dYe; + CCTK_REAL DEnt; vec dBvec; // Default constructor, no initialization. @@ -24,9 +25,24 @@ struct cons_vars { // Construct from single variables. CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline cons_vars( - CCTK_REAL dens_, vec mom_, CCTK_REAL tau_, CCTK_REAL dYe_, + CCTK_REAL dens_, vec mom_, CCTK_REAL tau_, CCTK_REAL dYe_, CCTK_REAL DEnt_, vec dBvec_) - : dens{dens_}, mom{mom_}, tau{tau_}, dYe{dYe_}, dBvec{dBvec_} {} + : dens{dens_}, mom{mom_}, tau{tau_}, dYe{dYe_}, DEnt{DEnt_}, dBvec{dBvec_} {} + + /// Copy assignment + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline cons_vars & + operator=(const cons_vars &other) { + if (this == &other) + return *this; // Handle self-assignment + // Copy data members from 'other' to 'this' + dens = other.dens; + mom = other.mom; + tau = other.tau; + dYe = other.dYe; + DEnt = other.DEnt; + dBvec = other.dBvec; + return *this; + } CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline void from_prim(const prim_vars &pv, const smat &g) { @@ -38,7 +54,8 @@ struct cons_vars { const vec &v_up = pv.vel; const vec v_low = calc_contraction(g, v_up); - const CCTK_REAL w_lorentz = calc_wlorentz(v_low, v_up); + const CCTK_REAL w_lorentz = pv.w_lor;; + //const CCTK_REAL w_lorentz = calc_wlorentz(v_low, v_up); /* Computing B_j */ const vec &B_up = pv.Bvec; @@ -68,11 +85,12 @@ struct cons_vars { dBvec = sqrt_detg * pv.Bvec; dYe = dens * pv.Ye; + DEnt = dens * pv.entropy; } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void scatter(CCTK_REAL &dens_, CCTK_REAL &momx_, CCTK_REAL &momy_, - CCTK_REAL &momz_, CCTK_REAL &tau_, CCTK_REAL &dYe_, + CCTK_REAL &momz_, CCTK_REAL &tau_, CCTK_REAL &dYe_, CCTK_REAL &DEnt_, CCTK_REAL &dBvecx_, CCTK_REAL &dBvecy_, CCTK_REAL &dBvecz_) const { dens_ = dens; @@ -81,13 +99,14 @@ struct cons_vars { momz_ = mom(2); tau_ = tau; dYe_ = dYe; + DEnt_ = DEnt; dBvecx_ = dBvec(0); dBvecy_ = dBvec(1); dBvecz_ = dBvec(2); } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void set_to_nan() { - dens = mom(0) = mom(1) = mom(2) = tau = dYe = dBvec(0) = dBvec(1) = + dens = mom(0) = mom(1) = mom(2) = tau = dYe = DEnt = dBvec(0) = dBvec(1) = dBvec(2) = std::numeric_limits::quiet_NaN(); } }; diff --git a/Con2PrimFactory/src/prims.hxx b/Con2PrimFactory/src/prims.hxx index c39a26c7..baaebbb3 100644 --- a/Con2PrimFactory/src/prims.hxx +++ b/Con2PrimFactory/src/prims.hxx @@ -16,6 +16,16 @@ struct prim_vars { CCTK_REAL eps; CCTK_REAL Ye; CCTK_REAL press; + CCTK_REAL entropy; +// entropy refers to the "evolved" entropy. +// However, note that this quantity is +// not necessarily the "physical" entropy. This depends on the EOS, e.g. +// for the ideal gas we have entropy = p rho^(-gamma). +// The distinction between "evolved" and "physical" entropy is made +// explicit in EOSX where functions with "entropy_..." and "kappa_..." +// refer to "physical" and "evolved" entropy, respectively. In this +// application thorn we always refer to entropy to describe +// the "evolved" entropy. vec vel; CCTK_REAL w_lor; vec Bvec; @@ -26,13 +36,31 @@ struct prim_vars { /// Construct from single variables. CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline prim_vars( - CCTK_REAL rho_, CCTK_REAL eps_, CCTK_REAL ye_, CCTK_REAL press_, + CCTK_REAL rho_, CCTK_REAL eps_, CCTK_REAL ye_, CCTK_REAL press_, CCTK_REAL entropy_, vec vel_, CCTK_REAL w_lor_, vec Bvec_) - : rho(rho_), eps(eps_), Ye(ye_), press(press_), vel(vel_), w_lor(w_lor_), + : rho(rho_), eps(eps_), Ye(ye_), press(press_), entropy(entropy_), vel(vel_), w_lor(w_lor_), Bvec(Bvec_){}; + /// Copy assignment + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline prim_vars & + operator=(const prim_vars &other) { + if (this == &other) + return *this; // Handle self-assignment + // Copy data members from 'other' to 'this' + rho = other.rho; + eps = other.eps; + Ye = other.Ye; + press = other.press; + entropy = other.entropy; + vel = other.vel; + w_lor = other.w_lor; + Bvec = other.Bvec; + E = other.E; + return *this; + } + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void - scatter(CCTK_REAL &rho_, CCTK_REAL &eps_, CCTK_REAL &ye_, CCTK_REAL &press_, + scatter(CCTK_REAL &rho_, CCTK_REAL &eps_, CCTK_REAL &ye_, CCTK_REAL &press_, CCTK_REAL &entropy_, CCTK_REAL &velx_, CCTK_REAL &vely_, CCTK_REAL &velz_, CCTK_REAL &w_lor_, CCTK_REAL &Bvec_x_, CCTK_REAL &Bvec_y_, CCTK_REAL &Bvec_z_, CCTK_REAL &E_x_, CCTK_REAL &E_y_, @@ -41,6 +69,7 @@ struct prim_vars { eps_ = eps; ye_ = Ye; press_ = press; + entropy_ = entropy; velx_ = vel(0); vely_ = vel(1); velz_ = vel(2); @@ -54,7 +83,7 @@ struct prim_vars { } CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void set_to_nan() { - rho = eps = Ye = press = vel(0) = vel(1) = vel(2) = w_lor = Bvec(0) = + rho = eps = Ye = press = entropy = vel(0) = vel(1) = vel(2) = w_lor = Bvec(0) = Bvec(1) = Bvec(2) = E(0) = E(1) = E(2) = std::numeric_limits::quiet_NaN(); } diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx index c80b592c..e1dc4bbd 100644 --- a/Con2PrimFactory/src/test.cxx +++ b/Con2PrimFactory/src/test.cxx @@ -9,6 +9,7 @@ #include "c2p.hxx" #include "c2p_1DPalenzuela.hxx" #include "c2p_2DNoble.hxx" +#include "c2p_1DEntropy.hxx" #include "c2p_utils.hxx" @@ -25,40 +26,84 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { eos::range rgeps(0, 100), rgrho(1e-13, 20), rgye(0.5, 0.5); const eos_idealgas eos_th(2.0, 938.985, rgeps, rgrho, rgye); + // Set atmo values + const CCTK_REAL rho_atmo = 1e-10; + const CCTK_REAL eps_atmo = 1e-8; + const CCTK_REAL Ye_atmo = 0.5; + const CCTK_REAL press_atmo = eos_th.press_from_valid_rho_eps_ye(rho_atmo,eps_atmo,Ye_atmo); + const CCTK_REAL entropy_atmo = eos_th.kappa_from_valid_rho_eps_ye(rho_atmo,eps_atmo,Ye_atmo); + // Setting up atmosphere - atmosphere atmo(1e-10, 1e-8, 0.5, 1e-8, 0.001); + const CCTK_REAL rho_atmo_cut = rho_atmo * (1 + 1.0e-3); + atmosphere atmo(rho_atmo, eps_atmo, Ye_atmo, press_atmo, entropy_atmo, rho_atmo_cut); // Metric const smat g{1.0, 0.0, 0.0, 1.0, 0.0, 1.0}; // xx, xy, xz, yy, yz, zz + // Set BH limiters + const CCTK_REAL alp_thresh = -1.; + const CCTK_REAL rho_BH = 1e20; + const CCTK_REAL eps_BH = 1e20; + const CCTK_REAL vwlim_BH = 1e20; + // Con2Prim objects - c2p_2DNoble c2p_Noble(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true); - c2p_1DPalenzuela c2p_Pal(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true); + // (eos_th, atmo, max_iter, c2p_tol + // alp_thresh, cons_error_limit, + // vw_lim, B_lim, rho_BH, eps_BH, vwlim_BH, + // Ye_lenient, use_z) + c2p_2DNoble c2p_Noble(eos_th, atmo, 100, 1e-8, alp_thresh, -1., 1, 1, rho_BH, eps_BH, vwlim_BH, true, false); + c2p_1DPalenzuela c2p_Pal(eos_th, atmo, 100, 1e-8, alp_thresh, -1., 1, 1, rho_BH, eps_BH, vwlim_BH, true, false); + c2p_1DEntropy c2p_Ent(eos_th, atmo, 100, 1e-8, alp_thresh, -1., 1, 1, rho_BH, eps_BH, vwlim_BH, true, false); // Construct error report object: c2p_report rep_Noble; c2p_report rep_Pal; + c2p_report rep_Ent; + + // Set primitive seeds + const CCTK_REAL rho_in = 0.125; + const CCTK_REAL eps_in = 0.8; + const CCTK_REAL Ye_in = 0.5; + const CCTK_REAL press_in = eos_th.press_from_valid_rho_eps_ye(rho_in,eps_in,Ye_in); + const CCTK_REAL entropy_in = eos_th.kappa_from_valid_rho_eps_ye(rho_in,eps_in,Ye_in); + const vec vup_in = {0.0,0.0,0.0}; + const vec Bup_in = {0.5,-1.0,0.0}; + const vec vdown_in = calc_contraction(g,vup_in); + const CCTK_REAL wlor_in = calc_wlorentz(vdown_in,vup_in); prim_vars pv; - // rho(p.I), eps(p.I), dummy_Ye, press(p.I),v_up, wlor, Bup - prim_vars pv_seeds{0.125, 0.8, 0.5, 0.1, {0, 0, 0}, 1, {0.5, -1.0, 0}}; + // rho(p.I), eps(p.I), dummy_Ye, press(p.I), entropy, v_up, wlor, Bup + prim_vars pv_seeds{rho_in, eps_in, Ye_in, press_in, entropy_in, vup_in, wlor_in, Bup_in}; // cons_vars cv{dens(p.I), {momx(p.I), momy(p.I), momz(p.I)}, tau(p.I), - // dummy_dYe, {dBx(p.I), dBy(p.I), dBz(p.I)}}; + // dummy_dYe, DEnt, {dBx(p.I), dBy(p.I), dBz(p.I)}}; + + cons_vars cv_Noble; + cons_vars cv_Pal; + cons_vars cv_Ent; + cons_vars cv_all; - cons_vars cv; - cv.from_prim(pv_seeds, g); + // C2P solve may modify cons_vars input + // Use the following to test C2Ps independently + cv_Noble.from_prim(pv_seeds, g); + cv_Pal.from_prim(pv_seeds, g); + cv_Ent.from_prim(pv_seeds, g); + + // Use the following to test C2Ps in sequence + // (As done in AsterX, the evolution thorn) + cv_all.from_prim(pv_seeds, g); // Testing C2P Noble CCTK_VINFO("Testing C2P Noble..."); - c2p_Noble.solve(eos_th, pv, pv_seeds, cv, g, rep_Noble); + c2p_Noble.solve(eos_th, pv, pv_seeds, cv_all, g, rep_Noble); printf("pv_seeds, pv: \n" "rho: %f, %f \n" "eps: %f, %f \n" "Ye: %f, %f \n" "press: %f, %f \n" + "entropy: %f, %f \n" "velx: %f, %f \n" "vely: %f, %f \n" "velz: %f, %f \n" @@ -66,7 +111,8 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { "By: %f, %f \n" "Bz: %f, %f \n", pv_seeds.rho, pv.rho, pv_seeds.eps, pv.eps, pv_seeds.Ye, pv.Ye, - pv_seeds.press, pv.press, pv_seeds.vel(0), pv.vel(0), pv_seeds.vel(1), + pv_seeds.press, pv.press, pv_seeds.entropy, pv.entropy, + pv_seeds.vel(0), pv.vel(0), pv_seeds.vel(1), pv.vel(1), pv_seeds.vel(2), pv.vel(2), pv_seeds.Bvec(0), pv.Bvec(0), pv_seeds.Bvec(1), pv.Bvec(1), pv_seeds.Bvec(2), pv.Bvec(2)); printf("cv: \n" @@ -78,9 +124,10 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { "dYe: %f \n" "dBx: %f \n" "dBy: %f \n" - "dBz: %f \n", - cv.dens, cv.tau, cv.mom(0), cv.mom(1), cv.mom(2), cv.dYe, cv.dBvec(0), - cv.dBvec(1), cv.dBvec(2)); + "dBz: %f \n" + "DEnt: %f \n", + cv_all.dens, cv_all.tau, cv_all.mom(0), cv_all.mom(1), cv_all.mom(2), cv_all.dYe, cv_all.dBvec(0), + cv_all.dBvec(1), cv_all.dBvec(2), cv_all.DEnt); /* assert(pv.rho == pv_seeds.rho); assert(pv.eps == pv_seeds.eps); @@ -93,13 +140,15 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { // Testing C2P Palenzuela CCTK_VINFO("Testing C2P Palenzuela..."); - c2p_Pal.solve(eos_th, pv, pv_seeds, cv, g, rep_Pal); + //c2p_Pal.solve(eos_th, pv, pv_seeds, cv, g, rep_Pal); + c2p_Pal.solve(eos_th, pv, cv_all, g, rep_Pal); printf("pv_seeds, pv: \n" "rho: %f, %f \n" "eps: %f, %f \n" "Ye: %f, %f \n" "press: %f, %f \n" + "entropy: %f, %f \n" "velx: %f, %f \n" "vely: %f, %f \n" "velz: %f, %f \n" @@ -107,7 +156,8 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { "By: %f, %f \n" "Bz: %f, %f \n", pv_seeds.rho, pv.rho, pv_seeds.eps, pv.eps, pv_seeds.Ye, pv.Ye, - pv_seeds.press, pv.press, pv_seeds.vel(0), pv.vel(0), pv_seeds.vel(1), + pv_seeds.press, pv.press, pv_seeds.entropy, pv.entropy, + pv_seeds.vel(0), pv.vel(0), pv_seeds.vel(1), pv.vel(1), pv_seeds.vel(2), pv.vel(2), pv_seeds.Bvec(0), pv.Bvec(0), pv_seeds.Bvec(1), pv.Bvec(1), pv_seeds.Bvec(2), pv.Bvec(2)); printf("cv: \n" @@ -119,9 +169,10 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { "dYe: %f \n" "dBx: %f \n" "dBy: %f \n" - "dBz: %f \n", - cv.dens, cv.tau, cv.mom(0), cv.mom(1), cv.mom(2), cv.dYe, cv.dBvec(0), - cv.dBvec(1), cv.dBvec(2)); + "dBz: %f \n" + "DEnt: %f \n", + cv_all.dens, cv_all.tau, cv_all.mom(0), cv_all.mom(1), cv_all.mom(2), cv_all.dYe, cv_all.dBvec(0), + cv_all.dBvec(1), cv_all.dBvec(2), cv_all.DEnt); /* assert(pv.rho == pv_seeds.rho); assert(pv.eps == pv_seeds.eps); @@ -130,6 +181,49 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { assert(pv.Bvec == pv_seeds.Bvec); */ rep_Pal.debug_message(); + + // Testing C2P Entropy + CCTK_VINFO("Testing C2P Entropy..."); + c2p_Ent.solve(eos_th, pv, cv_all, g, rep_Ent); + + printf("pv_seeds, pv: \n" + "rho: %f, %f \n" + "eps: %f, %f \n" + "Ye: %f, %f \n" + "press: %f, %f \n" + "entropy: %f, %f \n" + "velx: %f, %f \n" + "vely: %f, %f \n" + "velz: %f, %f \n" + "Bx: %f, %f \n" + "By: %f, %f \n" + "Bz: %f, %f \n", + pv_seeds.rho, pv.rho, pv_seeds.eps, pv.eps, pv_seeds.Ye, pv.Ye, + pv_seeds.press, pv.press, pv_seeds.entropy, pv.entropy, + pv_seeds.vel(0), pv.vel(0), pv_seeds.vel(1), + pv.vel(1), pv_seeds.vel(2), pv.vel(2), pv_seeds.Bvec(0), pv.Bvec(0), + pv_seeds.Bvec(1), pv.Bvec(1), pv_seeds.Bvec(2), pv.Bvec(2)); + printf("cv: \n" + "dens: %f \n" + "tau: %f \n" + "momx: %f \n" + "momy: %f \n" + "momz: %f \n" + "dYe: %f \n" + "dBx: %f \n" + "dBy: %f \n" + "dBz: %f \n" + "DEnt: %f \n", + cv_all.dens, cv_all.tau, cv_all.mom(0), cv_all.mom(1), cv_all.mom(2), cv_all.dYe, cv_all.dBvec(0), + cv_all.dBvec(1), cv_all.dBvec(2), cv_all.DEnt); + /* + assert(pv.rho == pv_seeds.rho); + assert(pv.eps == pv_seeds.eps); + assert(pv.press == pv_seeds.press); + assert(pv.vel == pv_seeds.vel); + assert(pv.Bvec == pv_seeds.Bvec); + */ + rep_Ent.debug_message(); } } // namespace Con2PrimFactory diff --git a/EOSX/.README.swp b/EOSX/.README.swp deleted file mode 100644 index fcc2e49a..00000000 Binary files a/EOSX/.README.swp and /dev/null differ diff --git a/EOSX/param.ccl b/EOSX/param.ccl index dcee5339..ee1e09e6 100644 --- a/EOSX/param.ccl +++ b/EOSX/param.ccl @@ -35,8 +35,8 @@ CCTK_REAL rho_max "Validity region: max density" STEERABLE=RECOVER CCTK_REAL rho_min "Validity region: min density" STEERABLE=RECOVER { - 0:* :: "" -} 0.0 + (0:* :: "" +} 1.0e-100 CCTK_REAL eps_max "Validity region: max internal specific energy" STEERABLE=RECOVER { diff --git a/EOSX/src/eos_idealgas.hxx b/EOSX/src/eos_idealgas.hxx index e355bcd6..7807cb97 100644 --- a/EOSX/src/eos_idealgas.hxx +++ b/EOSX/src/eos_idealgas.hxx @@ -84,6 +84,34 @@ public: const CCTK_REAL ye ///< Electron fraction \f$ Y_e \f$ ) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + press_from_valid_rho_kappa_ye( + const CCTK_REAL rho, + const CCTK_REAL kappa, // p/rho^gamma + const CCTK_REAL ye + ) const; + + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + eps_from_valid_rho_kappa_ye( + const CCTK_REAL rho, + const CCTK_REAL kappa, // p/rho^gamma + const CCTK_REAL ye + ) const; + + // Note that kappa implements a generic thermodynamic quantity + // meant to describe the "evolved" entropy by an evolution/application + // thorn. + // The notion of the "evolved" entropy (kappa) might differ from the definition + // of the actual entropy (entropy_from..., see above) for different EOS, + // e.g. for the ideal gas EOS we have kappa = p * (rho)^(-gamma), + // where gamma is the adiabatic index. + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL + kappa_from_valid_rho_eps_ye( + const CCTK_REAL rho, + const CCTK_REAL eps, + const CCTK_REAL ye + ) const; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline range range_eps_from_valid_rho_ye( const CCTK_REAL rho, ///< Rest mass density \f$ \rho \f$ @@ -173,6 +201,40 @@ eos_idealgas::eps_from_valid_rho_temp_ye(const CCTK_REAL rho, return temp / temp_over_eps; } +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +eos_idealgas::press_from_valid_rho_kappa_ye( + const CCTK_REAL rho, + const CCTK_REAL kappa, // p/rho^gamma + const CCTK_REAL ye +) const { + return kappa*pow(rho,gamma); +} + +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +eos_idealgas::eps_from_valid_rho_kappa_ye( + const CCTK_REAL rho, + const CCTK_REAL kappa, // p/rho^gamma + const CCTK_REAL ye +) const { + return kappa*pow(rho,gamma-1.0)/(gamma-1.0); +}; + +// Note that kappa implements a generic thermodynamic quantity +// meant to describe the "evolved" entropy by an evolution/application +// thorn. +// The notion of the "evolved" entropy (kappa) might differ from the definition +// of the actual entropy (entropy_from..., see above) for different EOS, +// e.g. for the ideal gas EOS we have kappa = p * (rho)^(-gamma), +// where gamma is the adiabatic index. +CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL +eos_idealgas::kappa_from_valid_rho_eps_ye( + const CCTK_REAL rho, + const CCTK_REAL eps, + const CCTK_REAL ye +) const { + return (gamma-1.0)*eps*pow(rho,1.0-gamma); +}; + CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline eos::range eos_idealgas::range_eps_from_valid_rho_ye(const CCTK_REAL rho, const CCTK_REAL ye) const { diff --git a/FishboneMoncriefIDX/interface.ccl b/FishboneMoncriefIDX/interface.ccl index b2fb062a..d0edbe76 100644 --- a/FishboneMoncriefIDX/interface.ccl +++ b/FishboneMoncriefIDX/interface.ccl @@ -1,6 +1,6 @@ IMPLEMENTS: FishboneMoncriefIDX -INHERITS: ADMBaseX CarpetX HydroBaseX AsterX +INHERITS: ADMBaseX HydroBaseX AsterX #USES INCLUDE HEADER: fixmath.hxx USES INCLUDE HEADER: loop_device.hxx