From 3a0d48a9117657ba98adca5110c7838f5f56a4eb Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Mon, 15 May 2023 10:51:37 -0600 Subject: [PATCH 01/51] start pin based off of bondi --- inputs/standing_accretion_shock.pin | 112 ++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 inputs/standing_accretion_shock.pin diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin new file mode 100644 index 000000000..a8dcb4d0a --- /dev/null +++ b/inputs/standing_accretion_shock.pin @@ -0,0 +1,112 @@ +# © 2021. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National +# Nuclear Security Administration. All rights in the program are +# reserved by Triad National Security, LLC, and the U.S. Department of +# Energy/National Nuclear Security Administration. The Government is +# granted for itself and others acting on its behalf a nonexclusive, +# paid-up, irrevocable worldwide license in this material to reproduce, +# prepare derivative works, distribute copies to the public, perform +# publicly and display publicly, and to permit others to do so. + + +problem = standing_accretion_shock +ix1_bc = outflow +ox1_bc = outflow +ix2_bc = reflect +ox2_bc = reflect + + +problem_id = standing_accretion_shock # problem ID: basename of output filenames + + +variables = p.density, & + c.density, & + p.velocity, & + c.momentum, & + p.energy, & + c.energy, & + pressure, & + cs, & + p.ye, & + g.c.coord, & + g.n.coord + +file_type = hdf5 # Tabular data dump +dt = 10.0 # time increment between outputs + + +nlim = -1 # cycle limit +tlim = 100.0 # time limit +integrator = rk2 # time integration algorithm +ncycle_out = 1 # interval for stdout summary info +dt_init_fact = 1.e-6 + + +nghost = 4 +#refinement = adaptive +#numlevel = 3 + +nx1 = 256 # Number of zones in X1-direction +x1min = 0.59 # minimum value of X1 +x1max = 5 # maximum value of X1 +ix1_bc = user # Inner-X1 boundary condition flag +ox1_bc = user # Outer-X1 boundary condition flag + +nx2 = 1 # Number of zones in X2-direction +x2min = 0 # minimum value of X2 +x2max = 1 #3.14159265359 #1 # maximum value of X2. Pi +ix2_bc = user # Inner-X2 boundary condition flag +ox2_bc = user # Outer-X2 boundary condition flag + +nx3 = 1 # Number of zones in X3-direction +x3min = 0 # minimum value of X3 +x3max = 6.28318530718 # maximum value of X3. 2*pi +ix3_bc = periodic # Inner-X3 boundary condition flag +ox3_bc = periodic # Outer-X3 boundary condition flfgag + +num_threads = 1 # maximum number of OMP threads + +# +#nx1 = 512 +#nx2 = 1 +#nx3 = 1 + + +field = c.c.bulk.rho +method = derivative_order_1 +max_level = 3 + + +type = IdealGas +Gamma = 1.4 +Cv = 2.5 + + +hydro = true +he = false +3t = false +rad = false + + +xorder = 2 +cfl = 0.8 +riemann = llf +recon = linear +c2p_method = robust +c2p_max_iter = 100 +c2p_tol = 1.e-8 +Ye = false +mhd = false + + +a = 0 + + +derefine_poles = false + + +eps_dissociation = 0.3 # amount of specific “thermal energy” available after nuclear binding energy has been removed from eps for dissociated nuclei. +p1 = 1 # initial postshock pressure +rho1 = 1 # initial postshock density From f7b9b66e0f838295a6791fb75e0cafaad843a00a Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 12:49:08 -0600 Subject: [PATCH 02/51] update input parameters --- inputs/standing_accretion_shock.pin | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index a8dcb4d0a..8f02cc57c 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -107,6 +107,7 @@ a = 0 derefine_poles = false +Mdot = 1.e-2 # Msun / sc +rShock = 200 # km +rPNS = 1.3 # Msun eps_dissociation = 0.3 # amount of specific “thermal energy” available after nuclear binding energy has been removed from eps for dissociated nuclei. -p1 = 1 # initial postshock pressure -rho1 = 1 # initial postshock density From e4932186951a309c4bed04fb0c6dd956b053ea9a Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 13:16:23 -0600 Subject: [PATCH 03/51] template pgen from bondi --- inputs/standing_accretion_shock.pin | 2 +- src/pgen/standing_accretion_shock.cpp | 183 ++++++++++++++++++++++++++ 2 files changed, 184 insertions(+), 1 deletion(-) create mode 100644 src/pgen/standing_accretion_shock.cpp diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 8f02cc57c..68dbdc17d 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -110,4 +110,4 @@ derefine_poles = false Mdot = 1.e-2 # Msun / sc rShock = 200 # km rPNS = 1.3 # Msun -eps_dissociation = 0.3 # amount of specific “thermal energy” available after nuclear binding energy has been removed from eps for dissociated nuclei. +eps_ff_ke = 0.3 # fraction of the free-fall kinetic energy at the initial location of the shock diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp new file mode 100644 index 000000000..f0d810df1 --- /dev/null +++ b/src/pgen/standing_accretion_shock.cpp @@ -0,0 +1,183 @@ +// © 2021. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract +// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +// is operated by Triad National Security, LLC for the U.S. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#include "geometry/boyer_lindquist.hpp" +#include "geometry/mckinney_gammie_ryan.hpp" +#include "pgen/pgen.hpp" +#include "utils/error_checking.hpp" + +// namespace phoebus { + +namespace standing_accretion_shock { + +KOKKOS_INLINE_FUNCTION +void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { + using namespace Geometry; + Real trans[NDFULL][NDFULL]; + LinearAlgebra::SetZero(trans, NDFULL, NDFULL); + const Real idenom = 1.0 / (r * r - 2.0 * r + a * a); + trans[0][0] = 1.0; + trans[0][1] = 2.0 * r * idenom; + trans[1][1] = 1.0; + trans[2][2] = 1.0; + trans[3][1] = a * idenom; + trans[3][3] = 1.0; + LinearAlgebra::SetZero(ucon_ks, NDFULL); + SPACETIMELOOP2(mu, nu) { ucon_ks[mu] += trans[mu][nu] * ucon_bl[nu]; } +} + +void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + + PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), + "Problem \"standing_accretion_shock\" requires \"FMKS\" geometry!"); + + auto rc = pmb->meshblock_data.Get().get(); + + PackIndexMap imap; + auto v = rc->PackVariables( + {fluid_prim::density, fluid_prim::velocity, fluid_prim::energy, fluid_prim::bfield, + fluid_prim::ye, fluid_prim::pressure, fluid_prim::temperature, fluid_prim::gamma1}, + imap); + + const int irho = imap[fluid_prim::density].first; + const int ivlo = imap[fluid_prim::velocity].first; + const int ivhi = imap[fluid_prim::velocity].second; + const int ieng = imap[fluid_prim::energy].first; + const int ib_lo = imap[fluid_prim::bfield].first; + const int ib_hi = imap[fluid_prim::bfield].second; + const int iye = imap[fluid_prim::ye].second; + const int iprs = imap[fluid_prim::pressure].first; + const int itmp = imap[fluid_prim::temperature].first; + const int igm1 = imap[fluid_prim::gamma1].first; + + // this only works with ideal gases + const std::string eos_type = pin->GetString("eos", "type"); + PARTHENON_REQUIRE_THROWS(eos_type == "IdealGas", + "Standing Accretion Shock setup only works with ideal gas"); + const Real gam = pin->GetReal("eos", "Gamma"); + const Real Cv = pin->GetReal("eos", "Cv"); + const Real n = 1.0 / (gam - 1.0); + PARTHENON_REQUIRE_THROWS(std::fabs(n - Cv) < 1.e-12, "Bondi requires Cv = 1/(Gamma-1)"); + PARTHENON_REQUIRE_THROWS(std::fabs(gam - 1.4) < 1.e-12, "Bondi requires gamma = 1.4"); + const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 1.e-2); + const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); + const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 1.3); + const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); + + // Solution constants + const Real uc = std::sqrt(1.0 / (2. * rs)); + const Real Vc = -std::sqrt(std::pow(uc, 2) / (1. - 3. * std::pow(uc, 2))); + const Real Tc = -n * std::pow(Vc, 2) / ((n + 1.) * (n * std::pow(Vc, 2) - 1.)); + const Real C1 = uc * std::pow(rs, 2) * std::pow(Tc, n); + const Real C2 = + std::pow(1. + (1. + n) * Tc, 2) * + (1. - 2. / rs + std::pow(C1, 2) / (std::pow(rs, 4) * std::pow(Tc, 2 * n))); + + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + + auto &coords = pmb->coords; + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + + const Real a = pin->GetReal("geometry", "a"); + auto bl = Geometry::BoyerLindquist(a); + + auto geom = Geometry::GetCoordinateSystem(rc); + + // set up transformation stuff + auto gpkg = pmb->packages.Get("geometry"); + bool derefine_poles = gpkg->Param("derefine_poles"); + Real h = gpkg->Param("h"); + Real xt = gpkg->Param("xt"); + Real alpha = gpkg->Param("alpha"); + Real x0 = gpkg->Param("x0"); + Real smooth = gpkg->Param("smooth"); + auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + + pmb->par_for( + "Phoebus::ProblemGenerator::Bondi", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + Real x1 = coords.Xc<1>(k, j, i); + const Real x2 = coords.Xc<2>(k, j, i); + const Real x3 = coords.Xc<3>(k, j, i); + + Real r = tr.bl_radius(x1); + while (r < Rhor) { + x1 += coords.Dxc<1>(i); + r = tr.bl_radius(x1); + } + + Real eos_lambda[2]; + if (iye > 0) { + v(iye, k, j, i) = 0.5; + eos_lambda[0] = v(iye, k, j, i); + } + + v(itmp, k, j, i) = get_bondi_temp(r, n, C1, C2, Tc, rs); + v(irho, k, j, i) = std::pow(v(itmp, k, j, i), n); + v(ieng, k, j, i) = v(irho, k, j, i) * v(itmp, k, j, i) / (gam - 1.0); + v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy( + v(irho, k, j, i), v(ieng, k, j, i) / v(irho, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); + Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; + ucon_bl[1] = -C1 / (std::pow(v(itmp, k, j, i), n) * std::pow(r, 2)); + + Real gcov[4][4]; + const Real th = tr.bl_theta(x1, x2); + bl.SpacetimeMetric(0.0, r, th, x3, gcov); + Real AA = gcov[0][0]; + Real BB = 2. * (gcov[0][1] * ucon_bl[1] + gcov[0][2] * ucon_bl[2] + + gcov[0][3] * ucon_bl[3]); + Real CC = 1. + gcov[1][1] * ucon_bl[1] * ucon_bl[1] + + gcov[2][2] * ucon_bl[2] * ucon_bl[2] + + gcov[3][3] * ucon_bl[3] * ucon_bl[3] + + 2. * (gcov[1][2] * ucon_bl[1] * ucon_bl[2] + + gcov[1][3] * ucon_bl[1] * ucon_bl[3] + + gcov[2][3] * ucon_bl[2] * ucon_bl[3]); + Real discr = BB * BB - 4. * AA * CC; + ucon_bl[0] = (-BB - std::sqrt(discr)) / (2. * AA); + const Real W_bl = ucon_bl[0] * bl.Lapse(0.0, r, th, x3); + + Real ucon[4]; + tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); + + // ucon won't be properly normalized here if x1 is not consistent with i + // so renormalize + geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); + AA = gcov[0][0]; + BB = 2. * (gcov[0][1] * ucon[1] + gcov[0][2] * ucon[2] + gcov[0][3] * ucon[3]); + CC = 1. + gcov[1][1] * ucon[1] * ucon[1] + gcov[2][2] * ucon[2] * ucon[2] + + gcov[3][3] * ucon[3] * ucon[3] + + 2. * (gcov[1][2] * ucon[1] * ucon[2] + gcov[1][3] * ucon[1] * ucon[3] + + gcov[2][3] * ucon[2] * ucon[3]); + discr = BB * BB - 4. * AA * CC; + PARTHENON_REQUIRE(discr > 0, "discr < 0"); + ucon[0] = (-BB - std::sqrt(discr)) / (2. * AA); + + // now get three velocity + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real beta[3]; + geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); + Real W = lapse * ucon[0]; + for (int d = 0; d < 3; d++) { + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + } + }); + + fluid::PrimitiveToConserved(rc); +} + +} // namespace standing_accretion_shock From 19fa5227ca8d489436e6ce8f27a0764c8b013bca Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 13:39:51 -0600 Subject: [PATCH 04/51] convert to code units, update defaults --- inputs/standing_accretion_shock.pin | 4 ++-- src/pgen/standing_accretion_shock.cpp | 13 +++++++++++-- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 68dbdc17d..0c6216bb6 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -107,7 +107,7 @@ a = 0 derefine_poles = false -Mdot = 1.e-2 # Msun / sc +Mdot = 0.2 # Msun / sec rShock = 200 # km -rPNS = 1.3 # Msun +rPNS = 60 # km eps_ff_ke = 0.3 # fraction of the free-fall kinetic energy at the initial location of the shock diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index f0d810df1..84afdb6fc 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -15,6 +15,7 @@ #include "geometry/mckinney_gammie_ryan.hpp" #include "pgen/pgen.hpp" #include "utils/error_checking.hpp" +#include "phoebus_utils/unit_conversions.hpp" // namespace phoebus { @@ -69,9 +70,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real n = 1.0 / (gam - 1.0); PARTHENON_REQUIRE_THROWS(std::fabs(n - Cv) < 1.e-12, "Bondi requires Cv = 1/(Gamma-1)"); PARTHENON_REQUIRE_THROWS(std::fabs(gam - 1.4) < 1.e-12, "Bondi requires gamma = 1.4"); - const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 1.e-2); + const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 1.3); + const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 60); const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); // Solution constants @@ -93,6 +94,14 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); + auto &unit_conv = + pmb->packages.Get("phoebus")->Param("unit_conv"); + + // convert to CGS then to code units + Mdot *= ( (solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode() ); + rShock *= (1.e5 * unit_conv.GetLengthCGSToCode() ); + rPNS *= (1.e5 * unit_conv.GetLengthCGSToCode() ); + auto geom = Geometry::GetCoordinateSystem(rc); // set up transformation stuff From 2b951fcaa6b53cedac4e2b790507992004d9591a Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 13:48:28 -0600 Subject: [PATCH 05/51] propagate PROBLEM --- src/CMakeLists.txt | 1 + src/pgen/pgen.hpp | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d1e558d50..e420b42c9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -99,6 +99,7 @@ add_executable(phoebus pgen/shock_tube.cpp pgen/thin_cooling.cpp pgen/sedov.cpp + pgen/standing_accretion_shock.cpp pgen/torus.cpp pgen/tov.cpp diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 2df825235..616eaf2de 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -47,7 +47,8 @@ using namespace parthenon::package::prelude; PROBLEM(torus) \ PROBLEM(p2c2p) \ PROBLEM(tov) \ - PROBLEM(homologous) + PROBLEM(homologous) \ + PROBLEM(standing_accretion_shock) // if you need problem-specific modifications to inputs, add the name here #define FOREACH_MODIFIER \ From 436d66c7744638f9952b4d9741ebf609672c9f84 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 15:57:13 -0600 Subject: [PATCH 06/51] edit --- inputs/standing_accretion_shock.pin | 4 ++-- src/pgen/standing_accretion_shock.cpp | 11 +---------- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 0c6216bb6..a45b29cdf 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -49,8 +49,8 @@ nghost = 4 #numlevel = 3 nx1 = 256 # Number of zones in X1-direction -x1min = 0.59 # minimum value of X1 -x1max = 5 # maximum value of X1 +x1min = 0 # minimum value of X1 (overwritten by rPNS) +x1max = 300 # maximum value of X1 (in km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 84afdb6fc..5e893d8df 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -75,15 +75,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 60); const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); - // Solution constants - const Real uc = std::sqrt(1.0 / (2. * rs)); - const Real Vc = -std::sqrt(std::pow(uc, 2) / (1. - 3. * std::pow(uc, 2))); - const Real Tc = -n * std::pow(Vc, 2) / ((n + 1.) * (n * std::pow(Vc, 2) - 1.)); - const Real C1 = uc * std::pow(rs, 2) * std::pow(Tc, n); - const Real C2 = - std::pow(1. + (1. + n) * Tc, 2) * - (1. - 2. / rs + std::pow(C1, 2) / (std::pow(rs, 4) * std::pow(Tc, 2 * n))); - IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); @@ -115,7 +106,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); pmb->par_for( - "Phoebus::ProblemGenerator::Bondi", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real x1 = coords.Xc<1>(k, j, i); const Real x2 = coords.Xc<2>(k, j, i); From d1da9c5b41502871776b7ba05b7ae6ef295747c2 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 17 May 2023 16:04:23 -0600 Subject: [PATCH 07/51] formatting ig!! --- src/pgen/pgen.hpp | 2 +- src/pgen/standing_accretion_shock.cpp | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 616eaf2de..048edd0cf 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -48,7 +48,7 @@ using namespace parthenon::package::prelude; PROBLEM(p2c2p) \ PROBLEM(tov) \ PROBLEM(homologous) \ - PROBLEM(standing_accretion_shock) + PROBLEM(standing_accretion_shock) // if you need problem-specific modifications to inputs, add the name here #define FOREACH_MODIFIER \ diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 5e893d8df..1e7b3936c 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -14,8 +14,8 @@ #include "geometry/boyer_lindquist.hpp" #include "geometry/mckinney_gammie_ryan.hpp" #include "pgen/pgen.hpp" -#include "utils/error_checking.hpp" #include "phoebus_utils/unit_conversions.hpp" +#include "utils/error_checking.hpp" // namespace phoebus { @@ -73,7 +73,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 60); - const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); + const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -89,9 +89,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { pmb->packages.Get("phoebus")->Param("unit_conv"); // convert to CGS then to code units - Mdot *= ( (solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode() ); - rShock *= (1.e5 * unit_conv.GetLengthCGSToCode() ); - rPNS *= (1.e5 * unit_conv.GetLengthCGSToCode() ); + Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); + rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); + rPNS *= (1.e5 * unit_conv.GetLengthCGSToCode()); auto geom = Geometry::GetCoordinateSystem(rc); @@ -106,8 +106,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); pmb->par_for( - "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int k, const int j, const int i) { + "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, + ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { Real x1 = coords.Xc<1>(k, j, i); const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); From b3a68ff23a2dd99842a18974bbbc98bfe9656eae Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Mon, 22 May 2023 10:58:46 -0600 Subject: [PATCH 08/51] more edits --- src/pgen/standing_accretion_shock.cpp | 31 ++++++++++++++++----------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 1e7b3936c..8b33a1cdb 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -61,19 +61,17 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const int itmp = imap[fluid_prim::temperature].first; const int igm1 = imap[fluid_prim::gamma1].first; - // this only works with ideal gases const std::string eos_type = pin->GetString("eos", "type"); - PARTHENON_REQUIRE_THROWS(eos_type == "IdealGas", - "Standing Accretion Shock setup only works with ideal gas"); + const EosType eos_type_enum = (eos_type == "IdealGas") ? IdealGas : StellarCollapse; + PARTHENON_REQUIRE_THROWS(eos_type == "IdealGas" || eos_type == "StellarCollapse", + "Standing Accretion Shock setup only works with ideal gas or stellar collapse EOS"); const Real gam = pin->GetReal("eos", "Gamma"); const Real Cv = pin->GetReal("eos", "Cv"); - const Real n = 1.0 / (gam - 1.0); - PARTHENON_REQUIRE_THROWS(std::fabs(n - Cv) < 1.e-12, "Bondi requires Cv = 1/(Gamma-1)"); - PARTHENON_REQUIRE_THROWS(std::fabs(gam - 1.4) < 1.e-12, "Bondi requires gamma = 1.4"); + const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 60); - const Real eps_ff_ke = pin->GetOrAddReal("standing_accretion_shock", "eps_ff_ke", 0.3); + const Real eps0 = pin->GetOrAddReal("standing_accretion_shock", "eps0", 1.e-20); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -112,11 +110,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); - Real r = tr.bl_radius(x1); - while (r < Rhor) { - x1 += coords.Dxc<1>(i); - r = tr.bl_radius(x1); - } + Real r = tr.bl_radius(x1); // r = e^(x1min) + //const Real lapse0 = geom.Lapse(CellLocation::Cent, rShock, j, i); + + //if (r < rShock) { + // const Real W0 = 1. / lapse0 + // const Real rho0 = Mdot / (4. * PI * std::pow(rShock,2) * W0 * std::abs(vr0)) // code units + // v(irho, k, j, i) = rho0; + // v(ieng, k, j, i) = eps0; + // } Real eos_lambda[2]; if (iye > 0) { @@ -124,6 +126,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_lambda[0] = v(iye, k, j, i); } + // set initial quantities? v(itmp, k, j, i) = get_bondi_temp(r, n, C1, C2, Tc, rs); v(irho, k, j, i) = std::pow(v(itmp, k, j, i), n); v(ieng, k, j, i) = v(irho, k, j, i) * v(itmp, k, j, i) / (gam - 1.0); @@ -132,7 +135,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; + + // set initial radial velocity? + Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; ucon_bl[1] = -C1 / (std::pow(v(itmp, k, j, i), n) * std::pow(r, 2)); Real gcov[4][4]; From 721ba2d7cc0e2696fe371ba59b293afef444c4a9 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 23 May 2023 18:21:05 -0600 Subject: [PATCH 09/51] update pin defaults --- inputs/standing_accretion_shock.pin | 31 ++++++++++++++++------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index a45b29cdf..ae65a82be 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -22,23 +22,23 @@ problem_id = standing_accretion_shock # problem ID: basename of output fi variables = p.density, & - c.density, & + c.density, & p.velocity, & - c.momentum, & + c.momentum, & p.energy, & c.energy, & - pressure, & + pressure, & cs, & p.ye, & g.c.coord, & g.n.coord file_type = hdf5 # Tabular data dump -dt = 10.0 # time increment between outputs +dt = 1.e-3 # time increment between outputs nlim = -1 # cycle limit -tlim = 100.0 # time limit +tlim = 1.0 # time limit integrator = rk2 # time integration algorithm ncycle_out = 1 # interval for stdout summary info dt_init_fact = 1.e-6 @@ -49,14 +49,14 @@ nghost = 4 #numlevel = 3 nx1 = 256 # Number of zones in X1-direction -x1min = 0 # minimum value of X1 (overwritten by rPNS) -x1max = 300 # maximum value of X1 (in km) +x1min = 3.44199 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 +x1max = 5.5 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag nx2 = 1 # Number of zones in X2-direction x2min = 0 # minimum value of X2 -x2max = 1 #3.14159265359 #1 # maximum value of X2. Pi +x2max = 1 #3.14159265359 #1 # maximum value of X2. Pi ix2_bc = user # Inner-X2 boundary condition flag ox2_bc = user # Outer-X2 boundary condition flag @@ -79,9 +79,9 @@ method = derivative_order_1 max_level = 3 -type = IdealGas -Gamma = 1.4 -Cv = 2.5 +type = StellarCollapse +filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 +use_sp5 = false hydro = true @@ -97,7 +97,7 @@ recon = linear c2p_method = robust c2p_max_iter = 100 c2p_tol = 1.e-8 -Ye = false +Ye = true mhd = false @@ -106,8 +106,11 @@ a = 0 derefine_poles = false + +scale_free = false +geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km +fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g + Mdot = 0.2 # Msun / sec rShock = 200 # km -rPNS = 60 # km -eps_ff_ke = 0.3 # fraction of the free-fall kinetic energy at the initial location of the shock From dcf9cd75fcbc346226f95a6f557f1bc7fede57e0 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 24 May 2023 15:41:14 -0600 Subject: [PATCH 10/51] add function to get T from rho target mach --- inputs/standing_accretion_shock.pin | 1 + src/pgen/pgen.cpp | 24 +++++++ src/pgen/pgen.hpp | 5 ++ src/pgen/standing_accretion_shock.cpp | 90 +++++++++++++++++---------- 4 files changed, 87 insertions(+), 33 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index ae65a82be..5913479ef 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -114,3 +114,4 @@ fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km +target_mach = 100 # target Mach number for upstream flow in preshock region diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 67a90336a..177cd03b4 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -84,4 +84,28 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real return rho * eroot; } +KOKKOS_FUNCTION +Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, + const Real Tmin, const Real Tmax, const Real Ye, + const Real vr0) { + root_find::RootFind root; + const Real epsilon = 1.e-10; + Real Troot = root.regula_falsi( + [&](const Real T) { + Real lambda[2]; + lambda[0] = Ye; + Real P = eos.PressureFromDensityTemperature(rho, T, lambda); + Real eps = eos.EnergyFromDensityTemperature(rho, T, lambda); + Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); + Real u = rho * eps; // convert eps / V to specific internal energy + Real w = rho + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P + Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w + Real mach = vr0 / cs; // radial component of preshock velocity + Real mach_res = mach - target_mach; + return mach_res; + }, + Tmin, Tmax, epsilon * mach_res, std::max(Tmin, epsilon)); + return Troot; +} + } // namespace phoebus diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 048edd0cf..94f72bf65 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -121,6 +121,11 @@ KOKKOS_FUNCTION Real energy_from_rho_P(const Microphysics::EOS::EOS &eos, const Real rho, const Real P, const Real emin, const Real emax, const Real Ye = 0.0); +KOKKOS_FUNCTION +Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, + const Real target_mach, const Real Tmin, const Real Tmax, + const Real Ye = 0.5, const Real vr0); + } // namespace phoebus #endif diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 8b33a1cdb..c7660ce4b 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -62,16 +62,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const int igm1 = imap[fluid_prim::gamma1].first; const std::string eos_type = pin->GetString("eos", "type"); - const EosType eos_type_enum = (eos_type == "IdealGas") ? IdealGas : StellarCollapse; - PARTHENON_REQUIRE_THROWS(eos_type == "IdealGas" || eos_type == "StellarCollapse", - "Standing Accretion Shock setup only works with ideal gas or stellar collapse EOS"); - const Real gam = pin->GetReal("eos", "Gamma"); - const Real Cv = pin->GetReal("eos", "Cv"); - - const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); - const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - const Real rPNS = pin->GetOrAddReal("standing_accretion_shock", "rPNS", 60); - const Real eps0 = pin->GetOrAddReal("standing_accretion_shock", "eps0", 1.e-20); + PARTHENON_REQUIRE_THROWS( + eos_type == "StellarCollapse", + "Standing Accretion Shock setup only works with StellarCollapse EOS"); + + const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", + 0.2); // pin in units of (Msun / sec) + const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", + 200); // pin in units of (km) + const Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", + 100); // pin in units of (dimensionless) IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -89,7 +89,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // convert to CGS then to code units Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - rPNS *= (1.e5 * unit_conv.GetLengthCGSToCode()); + const Real Tmin = + 11604525006.1598 * unit_conv.GetTemperatureCGSToCode(); // 1 MeV => K => code units + const Real Tmax = 2901131251539.96 * + unit_conv.GetTemperatureCGSToCode(); // 250 MeV => K => code units auto geom = Geometry::GetCoordinateSystem(rc); @@ -110,36 +113,57 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); - Real r = tr.bl_radius(x1); // r = e^(x1min) - //const Real lapse0 = geom.Lapse(CellLocation::Cent, rShock, j, i); - - //if (r < rShock) { - // const Real W0 = 1. / lapse0 - // const Real rho0 = Mdot / (4. * PI * std::pow(rShock,2) * W0 * std::abs(vr0)) // code units - // v(irho, k, j, i) = rho0; - // v(ieng, k, j, i) = eps0; - // } + Real r = tr.bl_radius(x1); // r = e^(x1min) + Real vel_rad; + // set Ye everywhere Real eos_lambda[2]; if (iye > 0) { v(iye, k, j, i) = 0.5; eos_lambda[0] = v(iye, k, j, i); } - // set initial quantities? - v(itmp, k, j, i) = get_bondi_temp(r, n, C1, C2, Tc, rs); - v(irho, k, j, i) = std::pow(v(itmp, k, j, i), n); - v(ieng, k, j, i) = v(irho, k, j, i) * v(itmp, k, j, i) / (gam - 1.0); - v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy( - v(irho, k, j, i), v(ieng, k, j, i) / v(irho, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); - - // set initial radial velocity? - Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; - ucon_bl[1] = -C1 / (std::pow(v(itmp, k, j, i), n) * std::pow(r, 2)); + if (r < rShock) { + const Real lapse0 = geom.Lapse( + CellLocation::Cent, k, j, + rShock); // cold, free falling, lapse equal to value at shock radius + const Real W0 = 1. / lapse0; + const Real vr0 = abs((std::sqrt(W0 - 1)) / (std::sqrt(W0))); + vel_rad = const Real rho0 = + Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + v(irho, k, j, i) = rho0; + Real T0 = temperature_from_rho_mach( + &eos, const Real rho0, const Real target_mach, const Real Tmin, + const Real Tmax, const Real Ye = 0.5, const Real vr0); + v(itmp, k, j, i) = T0; + v(ieng, k, j, i) = + rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); + + } else { + const Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); + const Real W0 = 1. / lapse0; + const Real vr0 = abs((std::sqrt(W0 - 1)) / (std::sqrt(W0))); + const Real rho0 = Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real gamma = 1.4; + const Real psi = std::pow(2, lapse1) * + ((gamma - 1) / gamma * (W0 - 1) / + W0); // here we assume epsilonND = 0 and accounted for in EoS + const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; + const Real rho1 = rho0*W0*(vr0/vr1); + v(irho, k, j, i) = rho1; + Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real Ye = 0.5, const Real vr1); + v(itmp, k, j, i) = T1; + v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + } + Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; Real gcov[4][4]; const Real th = tr.bl_theta(x1, x2); bl.SpacetimeMetric(0.0, r, th, x3, gcov); From 8c2055adb3bb3ee00571364b7882a35dfd6c2cdc Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 14:05:13 -0600 Subject: [PATCH 11/51] pass Ye correctly --- src/pgen/standing_accretion_shock.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index c7660ce4b..55dd86e51 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -134,7 +134,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i) = rho0; Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, - const Real Tmax, const Real Ye = 0.5, const Real vr0); + const Real Tmax, eos_lambda[0], const Real vr0); v(itmp, k, j, i) = T0; v(ieng, k, j, i) = rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); @@ -156,7 +156,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; const Real rho1 = rho0*W0*(vr0/vr1); v(irho, k, j, i) = rho1; - Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real Ye = 0.5, const Real vr1); + Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, eos_lambda[0], const Real vr1); v(itmp, k, j, i) = T1; v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); From efe17cc892833cc6f7344c2f5fff42b383bb4f69 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 14:18:18 -0600 Subject: [PATCH 12/51] change order of call function --- src/pgen/pgen.cpp | 4 ++-- src/pgen/pgen.hpp | 2 +- src/pgen/standing_accretion_shock.cpp | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 177cd03b4..28cad69de 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -86,8 +86,8 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real KOKKOS_FUNCTION Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real Ye, - const Real vr0) { + const Real Tmin, const Real Tmax, const Real vr0, + const Real Ye) { root_find::RootFind root; const Real epsilon = 1.e-10; Real Troot = root.regula_falsi( diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 94f72bf65..243caa375 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -124,7 +124,7 @@ Real energy_from_rho_P(const Microphysics::EOS::EOS &eos, const Real rho, const KOKKOS_FUNCTION Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, const Real target_mach, const Real Tmin, const Real Tmax, - const Real Ye = 0.5, const Real vr0); + const Real vr0, const Real Ye = 0.5); } // namespace phoebus diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 55dd86e51..727a8bd49 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -128,13 +128,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { CellLocation::Cent, k, j, rShock); // cold, free falling, lapse equal to value at shock radius const Real W0 = 1. / lapse0; - const Real vr0 = abs((std::sqrt(W0 - 1)) / (std::sqrt(W0))); + const Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); vel_rad = const Real rho0 = Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); v(irho, k, j, i) = rho0; Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, - const Real Tmax, eos_lambda[0], const Real vr0); + const Real Tmax, const Real vr0, eos_lambda[0]); v(itmp, k, j, i) = T0; v(ieng, k, j, i) = rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); @@ -156,7 +156,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; const Real rho1 = rho0*W0*(vr0/vr1); v(irho, k, j, i) = rho1; - Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, eos_lambda[0], const Real vr1); + Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); v(itmp, k, j, i) = T1; v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); From 12e32440ff9ea69664b42c6d68d266549cf87da9 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 15:11:38 -0600 Subject: [PATCH 13/51] compute initial radial velocity assuming it has some r dependence that is yet TBD --- src/pgen/pgen.cpp | 5 +- src/pgen/standing_accretion_shock.cpp | 141 ++++++++++++++------------ 2 files changed, 77 insertions(+), 69 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 28cad69de..58ee978ca 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -101,10 +101,9 @@ Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target Real w = rho + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = vr0 / cs; // radial component of preshock velocity - Real mach_res = mach - target_mach; - return mach_res; + return mach - target_mach; }, - Tmin, Tmax, epsilon * mach_res, std::max(Tmin, epsilon)); + Tmin, Tmax, epsilon * mach, std::max(Tmin, epsilon)); return Troot; } diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 727a8bd49..08983c7e2 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -116,6 +116,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real r = tr.bl_radius(x1); // r = e^(x1min) Real vel_rad; + const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); + const Real W0 = 1. / lapse0; + const Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); + const Real rho0 = Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real gamma = 1.4; + // set Ye everywhere Real eos_lambda[2]; if (iye > 0) { @@ -124,17 +130,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } if (r < rShock) { - const Real lapse0 = geom.Lapse( - CellLocation::Cent, k, j, - rShock); // cold, free falling, lapse equal to value at shock radius - const Real W0 = 1. / lapse0; - const Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - vel_rad = const Real rho0 = - Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - v(irho, k, j, i) = rho0; Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); + v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T0; v(ieng, k, j, i) = rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); @@ -144,69 +143,79 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - } else { - const Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); - const Real W0 = 1. / lapse0; - const Real vr0 = abs((std::sqrt(W0 - 1)) / (std::sqrt(W0))); - const Real rho0 = Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - const Real gamma = 1.4; - const Real psi = std::pow(2, lapse1) * - ((gamma - 1) / gamma * (W0 - 1) / - W0); // here we assume epsilonND = 0 and accounted for in EoS - const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; - const Real rho1 = rho0*W0*(vr0/vr1); - v(irho, k, j, i) = rho1; - Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); - v(itmp, k, j, i) = T1; - v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - } + vel_rad = vr0; + Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; + Real gcov[4][4]; + const Real th = tr.bl_theta(x1, x2); + bl.SpacetimeMetric(0.0, r, th, x3, gcov); + ucon_bl[0] = ucon_norm(ucon_bl, gcov); + + Real ucon[4]; + tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); + geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); + ucon[0] = ucon_norm(ucon, gcov); - Real ucon_bl[] = {0.0, 0.0, 0.0, 0.0}; - Real gcov[4][4]; - const Real th = tr.bl_theta(x1, x2); - bl.SpacetimeMetric(0.0, r, th, x3, gcov); - Real AA = gcov[0][0]; - Real BB = 2. * (gcov[0][1] * ucon_bl[1] + gcov[0][2] * ucon_bl[2] + - gcov[0][3] * ucon_bl[3]); - Real CC = 1. + gcov[1][1] * ucon_bl[1] * ucon_bl[1] + - gcov[2][2] * ucon_bl[2] * ucon_bl[2] + - gcov[3][3] * ucon_bl[3] * ucon_bl[3] + - 2. * (gcov[1][2] * ucon_bl[1] * ucon_bl[2] + - gcov[1][3] * ucon_bl[1] * ucon_bl[3] + - gcov[2][3] * ucon_bl[2] * ucon_bl[3]); - Real discr = BB * BB - 4. * AA * CC; - ucon_bl[0] = (-BB - std::sqrt(discr)) / (2. * AA); - const Real W_bl = ucon_bl[0] * bl.Lapse(0.0, r, th, x3); - - Real ucon[4]; - tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); - - // ucon won't be properly normalized here if x1 is not consistent with i - // so renormalize - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); - AA = gcov[0][0]; - BB = 2. * (gcov[0][1] * ucon[1] + gcov[0][2] * ucon[2] + gcov[0][3] * ucon[3]); - CC = 1. + gcov[1][1] * ucon[1] * ucon[1] + gcov[2][2] * ucon[2] * ucon[2] + - gcov[3][3] * ucon[3] * ucon[3] + - 2. * (gcov[1][2] * ucon[1] * ucon[2] + gcov[1][3] * ucon[1] * ucon[3] + - gcov[2][3] * ucon[2] * ucon[3]); - discr = BB * BB - 4. * AA * CC; - PARTHENON_REQUIRE(discr > 0, "discr < 0"); - ucon[0] = (-BB - std::sqrt(discr)) / (2. * AA); - - // now get three velocity - const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); - Real beta[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real W = lapse * ucon[0]; - for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real beta[3]; + geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); + Real W = lapse * ucon[0]; + + // finally compute three-velocity + for (int d = 0; d < 3; d++) { + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + } + + } else { + Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); + Real psi = std::pow(2, lapse1) * ((gamma - 1) / gamma * (W0 - 1) / W0); + Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; + Real rho1 = rho0*W0*(vr0/vr1); + Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); + v(irho, k, j, i) = rho1; + v(itmp, k, j, i) = T1; + v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + + vel_rad = vr1; + Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; + Real gcov[4][4]; + const Real th = tr.bl_theta(x1, x2); + bl.SpacetimeMetric(0.0, r, th, x3, gcov); + ucon_bl[0] = ucon_norm(ucon_bl, gcov); + + Real ucon[4]; + tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); + geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); + ucon[0] = ucon_norm(ucon, gcov); + + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real beta[3]; + geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); + Real W = lapse * ucon[0]; + + // finally compute three-velocity + for (int d = 0; d < 3; d++) { + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + } } }); fluid::PrimitiveToConserved(rc); } +KOKKOS_FUNCTION +Real ucon_norm(Real ucon[4], Real gcov[4][4]) { + Real AA = gcov[0][0]; + Real BB = 2. * (gcov[0][1] * ucon[1] + gcov[0][2] * ucon[2] + gcov[0][3] * ucon[3]); + Real CC = 1. + gcov[1][1] * ucon[1] * ucon[1] + gcov[2][2] * ucon[2] * ucon[2] + + gcov[3][3] * ucon[3] * ucon[3] + + 2. * (gcov[1][2] * ucon[1] * ucon[2] + gcov[1][3] * ucon[1] * ucon[3] + + gcov[2][3] * ucon[2] * ucon[3]); + Real discr = BB * BB - 4. * AA * CC; + if (discr < 0) printf("discr = %g %g %g %g\n", discr, AA, BB, CC); + PARTHENON_REQUIRE(discr >= 0, "discr < 0"); + return (-BB - std::sqrt(discr)) / (2. * AA); +} + } // namespace standing_accretion_shock From 4c25e5f6072d9b517f48a9311b3d9f780c3c25a8 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 15:40:03 -0600 Subject: [PATCH 14/51] try with the function in the actual pgen file --- src/pgen/pgen.cpp | 22 ---------------------- src/pgen/pgen.hpp | 5 ----- src/pgen/standing_accretion_shock.cpp | 25 ++++++++++++++++++++++++- 3 files changed, 24 insertions(+), 28 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 58ee978ca..4f48bf3b3 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -83,28 +83,6 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real Real eroot = root.regula_falsi(res, emin, emax, 1.e-10 * P, emin - 1.e10); return rho * eroot; } - -KOKKOS_FUNCTION -Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, - const Real Ye) { - root_find::RootFind root; - const Real epsilon = 1.e-10; - Real Troot = root.regula_falsi( - [&](const Real T) { - Real lambda[2]; - lambda[0] = Ye; - Real P = eos.PressureFromDensityTemperature(rho, T, lambda); - Real eps = eos.EnergyFromDensityTemperature(rho, T, lambda); - Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); - Real u = rho * eps; // convert eps / V to specific internal energy - Real w = rho + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P - Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w - Real mach = vr0 / cs; // radial component of preshock velocity - return mach - target_mach; - }, - Tmin, Tmax, epsilon * mach, std::max(Tmin, epsilon)); - return Troot; } } // namespace phoebus diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 243caa375..048edd0cf 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -121,11 +121,6 @@ KOKKOS_FUNCTION Real energy_from_rho_P(const Microphysics::EOS::EOS &eos, const Real rho, const Real P, const Real emin, const Real emax, const Real Ye = 0.0); -KOKKOS_FUNCTION -Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, - const Real target_mach, const Real Tmin, const Real Tmax, - const Real vr0, const Real Ye = 0.5); - } // namespace phoebus #endif diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 08983c7e2..9ea407d94 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -119,7 +119,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); const Real W0 = 1. / lapse0; const Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - const Real rho0 = Mdot / (4. * PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); const Real gamma = 1.4; // set Ye everywhere @@ -130,6 +129,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } if (r < rShock) { + Real rho0 = Mdot / (4. * PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); @@ -218,4 +218,27 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]) { return (-BB - std::sqrt(discr)) / (2. * AA); } +KOKKOS_FUNCTION +Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, + const Real Tmin, const Real Tmax, const Real vr0, + const Real Ye) { + root_find::RootFind root; + const Real epsilon = 1.e-10; + Real Troot = root.regula_falsi( + [&](const Real T) { + Real lambda[2]; + lambda[0] = Ye; + Real P = eos.PressureFromDensityTemperature(rho, T, lambda); + Real eps = eos.EnergyFromDensityTemperature(rho, T, lambda); + Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); + Real u = rho * eps; // convert eps / V to specific internal energy + Real w = rho + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P + Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w + Real mach = vr0 / cs; // radial component of preshock velocity + return mach - target_mach; + }, + Tmin, Tmax, epsilon * mach, std::max(Tmin, epsilon)); + return Troot; +} + } // namespace standing_accretion_shock From a8bef382c5e33bc027340f783e6d0fa06fdc7414 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 15:53:42 -0600 Subject: [PATCH 15/51] extra bracket, make r dependence for lapse too --- src/pgen/pgen.cpp | 1 - src/pgen/standing_accretion_shock.cpp | 7 +++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 4f48bf3b3..67a90336a 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -83,6 +83,5 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real Real eroot = root.regula_falsi(res, emin, emax, 1.e-10 * P, emin - 1.e10); return rho * eroot; } -} } // namespace phoebus diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 9ea407d94..8c8601082 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -115,10 +115,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real r = tr.bl_radius(x1); // r = e^(x1min) Real vel_rad; - - const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real W0 = 1. / lapse0; - const Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); const Real gamma = 1.4; // set Ye everywhere @@ -129,6 +125,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } if (r < rShock) { + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); + Real W0 = 1. / lapse0; + Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, From c1bb862de6eafb8f5baf5e4f6ce5eb2f9deff178 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:03:39 -0600 Subject: [PATCH 16/51] format --- src/pgen/standing_accretion_shock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 8c8601082..5439cc645 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -125,7 +125,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } if (r < rShock) { - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * PI * std::pow(r, 2) * W0 * std::abs(vr0)); From a3a4aefd21e8b3b28296d547f41a83f614bb5366 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:17:34 -0600 Subject: [PATCH 17/51] forgot prototypes --- src/pgen/standing_accretion_shock.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 5439cc645..1bb68d5fc 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -37,6 +37,15 @@ void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { SPACETIMELOOP2(mu, nu) { ucon_ks[mu] += trans[mu][nu] * ucon_bl[nu]; } } +// Prototypes +// ---------------------------------------------------------------------- +KOKKOS_FUNCTION +Real ucon_norm(Real ucon[4], Real gcov[4][4]); +Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, + const Real Tmin, const Real Tmax, const Real vr0, + const Real Ye); +// ---------------------------------------------------------------------- + void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), From c5d44d2ad82ce1eff9fb0c6d500268104e1c573c Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:29:41 -0600 Subject: [PATCH 18/51] acutally pull in the eos lol, and pi = M_PI, TIL --- src/pgen/standing_accretion_shock.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 1bb68d5fc..d5e55746d 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -21,6 +21,9 @@ namespace standing_accretion_shock { +using Microphysics::EOS::EOS; +enum EosType { StellarCollapse }; + KOKKOS_INLINE_FUNCTION void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { using namespace Geometry; @@ -137,7 +140,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - Real rho0 = Mdot / (4. * PI * std::pow(r, 2) * W0 * std::abs(vr0)); + Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( &eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); From 53c11c9d84844551b1b0be4880c6c449b37b8360 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:37:59 -0600 Subject: [PATCH 19/51] pull in root find --- src/pgen/standing_accretion_shock.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index d5e55746d..839c09ee6 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -14,6 +14,7 @@ #include "geometry/boyer_lindquist.hpp" #include "geometry/mckinney_gammie_ryan.hpp" #include "pgen/pgen.hpp" +#include "phoebus_utils/root_find.hpp" #include "phoebus_utils/unit_conversions.hpp" #include "utils/error_checking.hpp" @@ -142,7 +143,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( - &eos, const Real rho0, const Real target_mach, const Real Tmin, + eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T0; @@ -181,7 +182,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real psi = std::pow(2, lapse1) * ((gamma - 1) / gamma * (W0 - 1) / W0); Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; Real rho1 = rho0*W0*(vr0/vr1); - Real T1 = temperature_from_rho_mach(&eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); + Real T1 = temperature_from_rho_mach(eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); v(irho, k, j, i) = rho1; v(itmp, k, j, i) = T1; v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); From 70df9c49d79e66eea3093231784166fa5aabcf39 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:42:54 -0600 Subject: [PATCH 20/51] correctly convert pin params --- src/pgen/standing_accretion_shock.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 839c09ee6..8405cde54 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -100,8 +100,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { pmb->packages.Get("phoebus")->Param("unit_conv"); // convert to CGS then to code units - Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); - rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); + Mdot_conv = + Mdot * ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); + rShock_conv = rShock * 1.e5 * unit_conv.GetLengthCGSToCode(); const Real Tmin = 11604525006.1598 * unit_conv.GetTemperatureCGSToCode(); // 1 MeV => K => code units const Real Tmax = 2901131251539.96 * @@ -137,11 +138,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_lambda[0] = v(iye, k, j, i); } - if (r < rShock) { + if (r < rShock_conv) { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + Real rho0 = Mdot_conv / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); From db4c8083ae1c5c603a28626b9defe96680e57617 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 16:48:27 -0600 Subject: [PATCH 21/51] remove const, use original conversions --- src/pgen/standing_accretion_shock.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 8405cde54..3a3481d4a 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -79,10 +79,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_type == "StellarCollapse", "Standing Accretion Shock setup only works with StellarCollapse EOS"); - const Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", - 0.2); // pin in units of (Msun / sec) - const Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", - 200); // pin in units of (km) + Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", + 0.2); // pin in units of (Msun / sec) + Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", + 200); // pin in units of (km) const Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); // pin in units of (dimensionless) @@ -100,9 +100,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { pmb->packages.Get("phoebus")->Param("unit_conv"); // convert to CGS then to code units - Mdot_conv = - Mdot * ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); - rShock_conv = rShock * 1.e5 * unit_conv.GetLengthCGSToCode(); + Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); + rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); const Real Tmin = 11604525006.1598 * unit_conv.GetTemperatureCGSToCode(); // 1 MeV => K => code units const Real Tmax = 2901131251539.96 * @@ -138,11 +137,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_lambda[0] = v(iye, k, j, i); } - if (r < rShock_conv) { + if (r < rShock) { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - Real rho0 = Mdot_conv / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( eos, const Real rho0, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); From d2e3aa3f7371dbd48866ea1203362b0043d8dfcc Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 25 May 2023 17:13:00 -0600 Subject: [PATCH 22/51] something still wrong with the call to singularity-eos idk --- src/pgen/standing_accretion_shock.cpp | 33 ++++++++++++++------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 3a3481d4a..3428cb109 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -45,9 +45,9 @@ void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { // ---------------------------------------------------------------------- KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); -Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, - const Real Ye); +Real temperature_from_rho_mach(const EOS &eos, const EosType eos_type, const Real rho, + const Real target_mach, const Real Tmin, const Real Tmax, + const Real vr0, const Real Ye); // ---------------------------------------------------------------------- void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { @@ -91,7 +91,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); auto &coords = pmb->coords; - auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + auto eos_h = pmb->packages.Get("eos")->Param("h.EOS"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); @@ -143,15 +144,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T0 = temperature_from_rho_mach( - eos, const Real rho0, const Real target_mach, const Real Tmin, - const Real Tmax, const Real vr0, eos_lambda[0]); + eos_h, eos_type_enum, const Real rho0, const Real target_mach, + const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T0; v(ieng, k, j, i) = - rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature( + rho0 * eos_h.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); + v(iprs, k, j, i) = eos_h.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(igm1, k, j, i) = eos_h.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); @@ -182,12 +183,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real psi = std::pow(2, lapse1) * ((gamma - 1) / gamma * (W0 - 1) / W0); Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; Real rho1 = rho0*W0*(vr0/vr1); - Real T1 = temperature_from_rho_mach(eos, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); + Real T1 = temperature_from_rho_mach(eos_h, eos_type_enum, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); v(irho, k, j, i) = rho1; v(itmp, k, j, i) = T1; - v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + v(ieng, k, j, i) = rho1 * eos_h.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); + v(iprs, k, j, i) = eos_h.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos_h.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); vel_rad = vr1; Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; @@ -231,9 +232,9 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]) { } KOKKOS_FUNCTION -Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, - const Real Ye) { +Real temperature_from_rho_mach(const EOS &eos, const EosType eos_type const Real rho, + const Real target_mach, const Real Tmin, const Real Tmax, + const Real vr0, const Real Ye) { root_find::RootFind root; const Real epsilon = 1.e-10; Real Troot = root.regula_falsi( From 5352f0a20e658ac9ca0fd8e738796f7263e5cd61 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 30 May 2023 15:42:22 -0600 Subject: [PATCH 23/51] try explicit res function instead, idk how to get lambda to work --- src/pgen/pgen.cpp | 38 ++++++++++++++ src/pgen/pgen.hpp | 4 ++ src/pgen/standing_accretion_shock.cpp | 72 ++++++++++----------------- 3 files changed, 68 insertions(+), 46 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 67a90336a..d2d740bc3 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -84,4 +84,42 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real return rho * eroot; } +class MachResidual { + public: + KOKKOS_INLINE_FUNCTION + MachResidual(const EOS &eos, const Real rho, const Real vr0, const Real target_mach, + const Real Ye) + : eos_(eos), rho_(rho), vr0_(vr0), target_mach_(target_mach) { + lambda_[0] = Ye; + } + KOKKOS_INLINE_FUNCTION + Real operator()(const Real T) { + Real P = eos_.PressureFromDensityTemperature(rho_, T, lambda_); + Real eps = eos_.InternalEnergyFromDensityTemperature(rho_, T, lambda_); + Real bmod = eos_.BulkModulusFromDensityTemperature(rho_, T, lambda_); + Real u = rho_ * eps; // convert eps / V to specific internal energy + Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P + Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w + Real mach = vr0_ / cs; // radial component of preshock velocity + return mach - target_mach_; + } + + private: + const EOS &eos_; + Real rho_, vr0_, target_mach_; + Real lambda_[2]; +}; + +KOKKOS_FUNCTION +Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, + const Real Tmin, const Real Tmax, const Real vr0, + const Real Ye) { + const Real epsilon = 1.e-10; + MachResidual res(eos, rho, vr0, target_mach, Ye); + root_find::RootFind root; + Real Troot = + root.regula_falsi(res, Tmin, Tmax, epsilon * target_mach, std::max(Tmin, epsilon)); + return Troot; +} + } // namespace phoebus diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 048edd0cf..0bbb314ce 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -121,6 +121,10 @@ KOKKOS_FUNCTION Real energy_from_rho_P(const Microphysics::EOS::EOS &eos, const Real rho, const Real P, const Real emin, const Real emax, const Real Ye = 0.0); +Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, + const Real target_mach, const Real Tmin, const Real Tmax, + const Real vr0, const Real Ye); + } // namespace phoebus #endif diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 3428cb109..2c3de133f 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -22,9 +22,6 @@ namespace standing_accretion_shock { -using Microphysics::EOS::EOS; -enum EosType { StellarCollapse }; - KOKKOS_INLINE_FUNCTION void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { using namespace Geometry; @@ -45,9 +42,6 @@ void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { // ---------------------------------------------------------------------- KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); -Real temperature_from_rho_mach(const EOS &eos, const EosType eos_type, const Real rho, - const Real target_mach, const Real Tmin, const Real Tmax, - const Real vr0, const Real Ye); // ---------------------------------------------------------------------- void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { @@ -91,8 +85,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); auto &coords = pmb->coords; - auto eos = pmb->packages.Get("eos")->Param("d.EOS"); - auto eos_h = pmb->packages.Get("eos")->Param("h.EOS"); + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); @@ -143,16 +136,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real W0 = 1. / lapse0; Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - Real T0 = temperature_from_rho_mach( - eos_h, eos_type_enum, const Real rho0, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, eos_lambda[0]); + Real T0 = phoebus::temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, + vr0, eos_lambda[0]); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T0; v(ieng, k, j, i) = - rho0 * eos_h.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); - v(iprs, k, j, i) = eos_h.PressureFromDensityTemperature( + rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos_h.BulkModulusFromDensityTemperature( + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); @@ -179,18 +171,29 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } } else { + // beyond rShock, quantities defined by value at rShock + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); + Real W0 = 1. / lapse0; + Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); + Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); Real psi = std::pow(2, lapse1) * ((gamma - 1) / gamma * (W0 - 1) / W0); - Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4.*psi)/2; - Real rho1 = rho0*W0*(vr0/vr1); - Real T1 = temperature_from_rho_mach(eos_h, eos_type_enum, const Real rho1, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr1, eos_lambda[0]); - v(irho, k, j, i) = rho1; + Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2; + Real rho1 = rho0 * W0 * (vr0 / vr1); + Real T1 = phoebus::temperature_from_rho_mach(eos, rho1, target_mach, Tmin, Tmax, + vr1, eos_lambda[0]); + v(irho, k, j, i) = rho1; v(itmp, k, j, i) = T1; - v(ieng, k, j, i) = rho1 * eos_h.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); - v(iprs, k, j, i) = eos_h.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos_h.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + v(ieng, k, j, i) = + rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); - vel_rad = vr1; + vel_rad = vr1; Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; Real gcov[4][4]; const Real th = tr.bl_theta(x1, x2); @@ -209,7 +212,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // finally compute three-velocity for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } } }); @@ -231,27 +234,4 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]) { return (-BB - std::sqrt(discr)) / (2. * AA); } -KOKKOS_FUNCTION -Real temperature_from_rho_mach(const EOS &eos, const EosType eos_type const Real rho, - const Real target_mach, const Real Tmin, const Real Tmax, - const Real vr0, const Real Ye) { - root_find::RootFind root; - const Real epsilon = 1.e-10; - Real Troot = root.regula_falsi( - [&](const Real T) { - Real lambda[2]; - lambda[0] = Ye; - Real P = eos.PressureFromDensityTemperature(rho, T, lambda); - Real eps = eos.EnergyFromDensityTemperature(rho, T, lambda); - Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); - Real u = rho * eps; // convert eps / V to specific internal energy - Real w = rho + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P - Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w - Real mach = vr0 / cs; // radial component of preshock velocity - return mach - target_mach; - }, - Tmin, Tmax, epsilon * mach, std::max(Tmin, epsilon)); - return Troot; -} - } // namespace standing_accretion_shock From 4ef165473a58c3f5fb9fcfd2324bd0e805d9febe Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 30 May 2023 16:20:45 -0600 Subject: [PATCH 24/51] try EOS generated Tmin,Tmax for bracket --- src/pgen/standing_accretion_shock.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 2c3de133f..abee9934b 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -86,6 +86,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &coords = pmb->coords; auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + auto Tmin = pmb->packages.Get("eos")->Param("T_min"); + auto Tmax = pmb->packages.Get("eos")->Param("T_max"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); @@ -96,10 +98,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // convert to CGS then to code units Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - const Real Tmin = - 11604525006.1598 * unit_conv.GetTemperatureCGSToCode(); // 1 MeV => K => code units - const Real Tmax = 2901131251539.96 * - unit_conv.GetTemperatureCGSToCode(); // 250 MeV => K => code units + // const Real Tmin = 1.16e9 * unit_conv.GetTemperatureCGSToCode(); // 0.1 MeV => K => + // code units const Real Tmax = 1.85e12 * unit_conv.GetTemperatureCGSToCode(); // 160 + // MeV => K => code units auto geom = Geometry::GetCoordinateSystem(rc); From 34e9c76faeb485e134c5d0b392e6d85b88d8b18f Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 31 May 2023 17:16:37 -0600 Subject: [PATCH 25/51] psi not being computed correctly. try this --- src/pgen/pgen.cpp | 7 +++-- src/pgen/standing_accretion_shock.cpp | 39 +++++++++++++++++++-------- 2 files changed, 33 insertions(+), 13 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index d2d740bc3..d4a9cc196 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -101,6 +101,8 @@ class MachResidual { Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = vr0_ / cs; // radial component of preshock velocity + // printf("Mach Residual Func produced: u, w, cs, Mach, vr0 = %g %g %g %g %g\n", u, w, + // cs, mach, vr0_); return mach - target_mach_; } @@ -114,11 +116,12 @@ KOKKOS_FUNCTION Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, const Real Ye) { - const Real epsilon = 1.e-10; + // printf("passing values of: rho, vr0, Mach_target, Ye = %g %g %g %g\n", rho, vr0, + // target_mach, Ye); MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; Real Troot = - root.regula_falsi(res, Tmin, Tmax, epsilon * target_mach, std::max(Tmin, epsilon)); + root.regula_falsi(res, Tmin, Tmax, 1.e-10 * target_mach, std::max(Tmin, 1.e-10)); return Troot; } diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index abee9934b..85851141a 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -88,6 +88,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); + printf("from SC: Tmin, Tmax = %g %g\n", Tmin, Tmax); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); @@ -98,9 +99,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // convert to CGS then to code units Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - // const Real Tmin = 1.16e9 * unit_conv.GetTemperatureCGSToCode(); // 0.1 MeV => K => - // code units const Real Tmax = 1.85e12 * unit_conv.GetTemperatureCGSToCode(); // 160 - // MeV => K => code units + Real Tmin_in = + 1.16e9 * unit_conv.GetTemperatureCGSToCode(); // 0.1 MeV => K => code units + Real Tmax_in = + 1.85e12 * unit_conv.GetTemperatureCGSToCode(); // 160 MeV => K => code units + printf("from input: Tmin_in, Tmax_in = %g %g\n", Tmin_in, Tmax_in); + const Real Mpns = 1.3 * solar_mass * unit_conv.GetMassCGSToCode(); auto geom = Geometry::GetCoordinateSystem(rc); @@ -123,7 +127,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real r = tr.bl_radius(x1); // r = e^(x1min) Real vel_rad; - const Real gamma = 1.4; + const Real gamma = 4. / 3.; // set Ye everywhere Real eos_lambda[2]; @@ -135,10 +139,20 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (r < rShock) { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; - Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); + Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; + // Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + Real alphasq = 1 - (2. * 1.3 / r); + Real psi = alphasq * ((gamma - 1) / gamma) * ((W0 - 1) / W0); + Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + Real rho1 = rho0 * W0 * (vr0 / vr1); + printf("inside shock: passing values of: r, rho0, vr0, vr1, Mach_target, " + "lapse0, psi = %g %g %g %g %g %g %g\n", + r, rho0, vr0, vr1, target_mach, lapse0, psi); Real T0 = phoebus::temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, - vr0, eos_lambda[0]); + vr1, eos_lambda[0]); + // printf("found T0 of: T0 at r, rho0, vr0, Mach_target = %g %g %g %g %g \n", + // T0, r, rho0, vr0, target_mach); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T0; v(ieng, k, j, i) = @@ -175,13 +189,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // beyond rShock, quantities defined by value at rShock Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; - Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); + Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; + // Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - - Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); - Real psi = std::pow(2, lapse1) * ((gamma - 1) / gamma * (W0 - 1) / W0); - Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2; + Real alphasq = 1 - (2. * Mpns / r); + Real psi = alphasq * ((gamma - 1) / gamma) * ((W0 - 1) / W0); + Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); + printf("outside shock: passing values of: r, rho0, vr0, vr1, Mach_target, " + "lapse0, psi = %g %g %g %g %g %g %g\n", + r, rho0, vr0, vr1, target_mach, lapse0, psi); Real T1 = phoebus::temperature_from_rho_mach(eos, rho1, target_mach, Tmin, Tmax, vr1, eos_lambda[0]); v(irho, k, j, i) = rho1; From 55ccd6994a7d9baa266de6c38d70ebf0087757fe Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Fri, 2 Jun 2023 11:23:02 -0600 Subject: [PATCH 26/51] correct pre/postshock assumptions, remove extra vars --- src/pgen/standing_accretion_shock.cpp | 75 +++++++++++---------------- 1 file changed, 31 insertions(+), 44 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 85851141a..2426e7652 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -85,10 +85,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); auto &coords = pmb->coords; - auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + auto eos = pmb->packages.Get("eos")->Param("h.EOS"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); - printf("from SC: Tmin, Tmax = %g %g\n", Tmin, Tmax); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); @@ -99,12 +98,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // convert to CGS then to code units Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - Real Tmin_in = - 1.16e9 * unit_conv.GetTemperatureCGSToCode(); // 0.1 MeV => K => code units - Real Tmax_in = - 1.85e12 * unit_conv.GetTemperatureCGSToCode(); // 160 MeV => K => code units - printf("from input: Tmin_in, Tmax_in = %g %g\n", Tmin_in, Tmax_in); - const Real Mpns = 1.3 * solar_mass * unit_conv.GetMassCGSToCode(); auto geom = Geometry::GetCoordinateSystem(rc); @@ -136,27 +129,21 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_lambda[0] = v(iye, k, j, i); } - if (r < rShock) { - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); + if (r > rShock) { + // preshock - 0 + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; - // Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - Real alphasq = 1 - (2. * 1.3 / r); - Real psi = alphasq * ((gamma - 1) / gamma) * ((W0 - 1) / W0); - Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - Real rho1 = rho0 * W0 * (vr0 / vr1); - printf("inside shock: passing values of: r, rho0, vr0, vr1, Mach_target, " - "lapse0, psi = %g %g %g %g %g %g %g\n", - r, rho0, vr0, vr1, target_mach, lapse0, psi); - Real T0 = phoebus::temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, - vr1, eos_lambda[0]); - // printf("found T0 of: T0 at r, rho0, vr0, Mach_target = %g %g %g %g %g \n", - // T0, r, rho0, vr0, target_mach); + Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + // printf("preshock (r>rShock): passing values of: r, rho0, vr0, lapse0, psi = + // %g %g %g %g\n",r, rho0, vr0, lapse0); + + Real T = phoebus::temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, + vr0, eos_lambda[0]); v(irho, k, j, i) = rho0; - v(itmp, k, j, i) = T0; + v(itmp, k, j, i) = T; v(ieng, k, j, i) = - rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T0, eos_lambda); + rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( @@ -175,36 +162,36 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); ucon[0] = ucon_norm(ucon, gcov); - const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + const Real lapsed = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real W = lapse * ucon[0]; + Real Wd = lapsed * ucon[0]; // finally compute three-velocity for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + v(ivlo + d, k, j, i) = ucon[d + 1] + Wd * beta[d] / lapsed; } } else { - // beyond rShock, quantities defined by value at rShock + // postshock - 1 Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; - // Real vr0 = abs((std::sqrt(W0 - 1.)) / (std::sqrt(W0))); - Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - Real alphasq = 1 - (2. * Mpns / r); - Real psi = alphasq * ((gamma - 1) / gamma) * ((W0 - 1) / W0); + Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + + Real alphasq = 1. - (2. / r); + Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); - printf("outside shock: passing values of: r, rho0, vr0, vr1, Mach_target, " - "lapse0, psi = %g %g %g %g %g %g %g\n", - r, rho0, vr0, vr1, target_mach, lapse0, psi); - Real T1 = phoebus::temperature_from_rho_mach(eos, rho1, target_mach, Tmin, Tmax, - vr1, eos_lambda[0]); + + const Real epsND = 0.003; // 2.7e16 ergs / g for M = 1.3Mpns + + // printf("postshock r Date: Fri, 2 Jun 2023 11:32:44 -0600 Subject: [PATCH 27/51] device not host --- src/pgen/standing_accretion_shock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 2426e7652..4c73540fb 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -85,7 +85,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); auto &coords = pmb->coords; - auto eos = pmb->packages.Get("eos")->Param("h.EOS"); + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); From 1e1b3be72716b3bac729877465c99696ac32cf07 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Fri, 2 Jun 2023 17:25:43 -0600 Subject: [PATCH 28/51] dont need transformation (i think), try to properly fill vel arrays --- src/pgen/standing_accretion_shock.cpp | 65 ++------------------------- 1 file changed, 3 insertions(+), 62 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 4c73540fb..0854830f9 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -22,28 +22,6 @@ namespace standing_accretion_shock { -KOKKOS_INLINE_FUNCTION -void bl_to_ks(const Real r, const Real a, Real *ucon_bl, Real *ucon_ks) { - using namespace Geometry; - Real trans[NDFULL][NDFULL]; - LinearAlgebra::SetZero(trans, NDFULL, NDFULL); - const Real idenom = 1.0 / (r * r - 2.0 * r + a * a); - trans[0][0] = 1.0; - trans[0][1] = 2.0 * r * idenom; - trans[1][1] = 1.0; - trans[2][2] = 1.0; - trans[3][1] = a * idenom; - trans[3][3] = 1.0; - LinearAlgebra::SetZero(ucon_ks, NDFULL); - SPACETIMELOOP2(mu, nu) { ucon_ks[mu] += trans[mu][nu] * ucon_bl[nu]; } -} - -// Prototypes -// ---------------------------------------------------------------------- -KOKKOS_FUNCTION -Real ucon_norm(Real ucon[4], Real gcov[4][4]); -// ---------------------------------------------------------------------- - void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), @@ -119,7 +97,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x3 = coords.Xc<3>(k, j, i); Real r = tr.bl_radius(x1); // r = e^(x1min) - Real vel_rad; const Real gamma = 4. / 3.; // set Ye everywhere @@ -150,18 +127,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - vel_rad = vr0; - Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; - Real gcov[4][4]; - const Real th = tr.bl_theta(x1, x2); - bl.SpacetimeMetric(0.0, r, th, x3, gcov); - ucon_bl[0] = ucon_norm(ucon_bl, gcov); - - Real ucon[4]; - tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); - ucon[0] = ucon_norm(ucon, gcov); - + Real ucon[4] = {0.0, vr0, 0.0, 0.0}; const Real lapsed = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); @@ -177,7 +143,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; - Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); Real alphasq = 1. - (2. / r); Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); @@ -198,18 +164,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - vel_rad = vr1; - Real ucon_bl[] = {0.0, vel_rad, 0.0, 0.0}; - Real gcov[4][4]; - const Real th = tr.bl_theta(x1, x2); - bl.SpacetimeMetric(0.0, r, th, x3, gcov); - ucon_bl[0] = ucon_norm(ucon_bl, gcov); - - Real ucon[4]; - tr.bl_to_fmks(x1, x2, x3, a, ucon_bl, ucon); - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); - ucon[0] = ucon_norm(ucon, gcov); - + Real ucon[4] = {0.0, vr1, 0.0, 0.0}; const Real lapsed = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); @@ -225,18 +180,4 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { fluid::PrimitiveToConserved(rc); } -KOKKOS_FUNCTION -Real ucon_norm(Real ucon[4], Real gcov[4][4]) { - Real AA = gcov[0][0]; - Real BB = 2. * (gcov[0][1] * ucon[1] + gcov[0][2] * ucon[2] + gcov[0][3] * ucon[3]); - Real CC = 1. + gcov[1][1] * ucon[1] * ucon[1] + gcov[2][2] * ucon[2] * ucon[2] + - gcov[3][3] * ucon[3] * ucon[3] + - 2. * (gcov[1][2] * ucon[1] * ucon[2] + gcov[1][3] * ucon[1] * ucon[3] + - gcov[2][3] * ucon[2] * ucon[3]); - Real discr = BB * BB - 4. * AA * CC; - if (discr < 0) printf("discr = %g %g %g %g\n", discr, AA, BB, CC); - PARTHENON_REQUIRE(discr >= 0, "discr < 0"); - return (-BB - std::sqrt(discr)) / (2. * AA); -} - } // namespace standing_accretion_shock From ba01a025adee863923fdb1335ea0cdf32adc148f Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 8 Jun 2023 12:35:27 -0600 Subject: [PATCH 29/51] move mach stuff into sasw pgen, add missing decorator, keep prints for now while debug root find failure --- src/pgen/pgen.cpp | 41 ------ src/pgen/pgen.hpp | 4 - src/pgen/standing_accretion_shock.cpp | 202 ++++++++++++++++++-------- 3 files changed, 143 insertions(+), 104 deletions(-) diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index d4a9cc196..67a90336a 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -84,45 +84,4 @@ Real energy_from_rho_P(const EOS &eos, const Real rho, const Real P, const Real return rho * eroot; } -class MachResidual { - public: - KOKKOS_INLINE_FUNCTION - MachResidual(const EOS &eos, const Real rho, const Real vr0, const Real target_mach, - const Real Ye) - : eos_(eos), rho_(rho), vr0_(vr0), target_mach_(target_mach) { - lambda_[0] = Ye; - } - KOKKOS_INLINE_FUNCTION - Real operator()(const Real T) { - Real P = eos_.PressureFromDensityTemperature(rho_, T, lambda_); - Real eps = eos_.InternalEnergyFromDensityTemperature(rho_, T, lambda_); - Real bmod = eos_.BulkModulusFromDensityTemperature(rho_, T, lambda_); - Real u = rho_ * eps; // convert eps / V to specific internal energy - Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P - Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w - Real mach = vr0_ / cs; // radial component of preshock velocity - // printf("Mach Residual Func produced: u, w, cs, Mach, vr0 = %g %g %g %g %g\n", u, w, - // cs, mach, vr0_); - return mach - target_mach_; - } - - private: - const EOS &eos_; - Real rho_, vr0_, target_mach_; - Real lambda_[2]; -}; - -KOKKOS_FUNCTION -Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, - const Real Ye) { - // printf("passing values of: rho, vr0, Mach_target, Ye = %g %g %g %g\n", rho, vr0, - // target_mach, Ye); - MachResidual res(eos, rho, vr0, target_mach, Ye); - root_find::RootFind root; - Real Troot = - root.regula_falsi(res, Tmin, Tmax, 1.e-10 * target_mach, std::max(Tmin, 1.e-10)); - return Troot; -} - } // namespace phoebus diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 0bbb314ce..048edd0cf 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -121,10 +121,6 @@ KOKKOS_FUNCTION Real energy_from_rho_P(const Microphysics::EOS::EOS &eos, const Real rho, const Real P, const Real emin, const Real emax, const Real Ye = 0.0); -Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, - const Real target_mach, const Real Tmin, const Real Tmax, - const Real vr0, const Real Ye); - } // namespace phoebus #endif diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 0854830f9..81f52ff46 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -12,20 +12,59 @@ // publicly, and to permit others to do so. #include "geometry/boyer_lindquist.hpp" -#include "geometry/mckinney_gammie_ryan.hpp" #include "pgen/pgen.hpp" #include "phoebus_utils/root_find.hpp" #include "phoebus_utils/unit_conversions.hpp" #include "utils/error_checking.hpp" -// namespace phoebus { - namespace standing_accretion_shock { +parthenon::constants::PhysicalConstants pc; +using Microphysics::EOS::EOS; + +class MachResidual { + public: + KOKKOS_FUNCTION + MachResidual(const EOS &eos, const Real rho, const Real vr0, const Real target_mach, + const Real Ye) + : eos_(eos), rho_(rho), vr0_(vr0), target_mach_(target_mach) { + lambda_[0] = Ye; + } + KOKKOS_INLINE_FUNCTION + Real operator()(const Real T) { + Real P = eos_.PressureFromDensityTemperature(rho_, T, lambda_); + Real eps = eos_.InternalEnergyFromDensityTemperature(rho_, T, lambda_); + Real bmod = eos_.BulkModulusFromDensityTemperature(rho_, T, lambda_); + Real u = rho_ * eps; // convert eps / V to specific internal energy + Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P + Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w + Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity + // printf("Mach Residual Func produced: u, w, cs, Mach, vr0 = %g %g %g %g %g\n", u, w, + // cs, mach, vr0_); + return mach - target_mach_; + } + + private: + const EOS &eos_; + Real rho_, vr0_, target_mach_; + Real lambda_[2]; +}; + +// protehypes---- +KOKKOS_FUNCTION +Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, + const Real target_mach, const Real Tmin, const Real Tmax, + const Real vr0, const Real Ye); +KOKKOS_FUNCTION +Real ucon_norm(Real ucon[4], Real gcov[4][4]); + +//-------------- + void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), - "Problem \"standing_accretion_shock\" requires \"FMKS\" geometry!"); + PARTHENON_REQUIRE( + typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::BoyerLindquist), + "Problem \"standing_accretion_shock\" requires \"BoyerLindquist\" geometry!"); auto rc = pmb->meshblock_data.Get().get(); @@ -51,12 +90,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_type == "StellarCollapse", "Standing Accretion Shock setup only works with StellarCollapse EOS"); - Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", - 0.2); // pin in units of (Msun / sec) - Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", - 200); // pin in units of (km) - const Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", - 100); // pin in units of (dimensionless) + Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); + Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); + const Real target_mach = + pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -73,21 +110,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &unit_conv = pmb->packages.Get("phoebus")->Param("unit_conv"); - // convert to CGS then to code units Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); + const Real MPNS = 1.3 * solar_mass * unit_conv.GetMassCGSToCode(); + const Real rs = + (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); - auto geom = Geometry::GetCoordinateSystem(rc); + printf("rShock (code units) = %g \n", rShock); - // set up transformation stuff - auto gpkg = pmb->packages.Get("geometry"); - bool derefine_poles = gpkg->Param("derefine_poles"); - Real h = gpkg->Param("h"); - Real xt = gpkg->Param("xt"); - Real alpha = gpkg->Param("alpha"); - Real x0 = gpkg->Param("x0"); - Real smooth = gpkg->Param("smooth"); - auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + auto geom = Geometry::GetCoordinateSystem(rc); pmb->par_for( "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, @@ -96,8 +127,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); - Real r = tr.bl_radius(x1); // r = e^(x1min) + Real r = std::abs(x1); const Real gamma = 4. / 3.; + const Real epsND = 0.003; // 2.7e16 ergs / g for M = 1.3Mpns + + // printf("0 (preshock region) x1, r, rho0, vr0, lapse0, W0 = %g %g %g %g %g + // %g\n",x1, r, rho0, vr0, lapse0, W0); // set Ye everywhere Real eos_lambda[2]; @@ -106,78 +141,127 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_lambda[0] = v(iye, k, j, i); } - if (r > rShock) { - // preshock - 0 - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); + // postshock - 1 + if (r < rShock) { + + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real W0 = 1. / lapse0; - Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; - Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - // printf("preshock (r>rShock): passing values of: r, rho0, vr0, lapse0, psi = - // %g %g %g %g\n",r, rho0, vr0, lapse0); + Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - Real T = phoebus::temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, - vr0, eos_lambda[0]); - v(irho, k, j, i) = rho0; - v(itmp, k, j, i) = T; - v(ieng, k, j, i) = - rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); + Real alphasq = 1. - (2. / r); + Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); + Real vr1 = -1. * (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + Real rho1 = rho0 * W0 * (vr0 / vr1); + + printf("inside shock (postshock region) r, rho1, vr1 = %g %g %g\n", r, rho1, + vr1); + // printf("inside shock (postshock region) r, alphasq, psi, rho0, vr0, rho1, + // vr1 = %g %g %g %g %g %g %g\n",r, alphasq, psi, rho0, vr0, rho1, vr1); + v(irho, k, j, i) = rho1; + v(ieng, k, j, i) = (W0 - 1. + epsND * (gamma - 1.)) / (gamma); + v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( + rho1, v(ieng, k, j, i) / rho1, eos_lambda); v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - Real ucon[4] = {0.0, vr0, 0.0, 0.0}; - const Real lapsed = geom.Lapse(CellLocation::Cent, k, j, i); + // construct contravariant 4-velocity + Real ucon[] = {0.0, vr1, 0.0, 0.0}; + + // initialize spacetime metric + Real gcov[4][4]; + geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); + ucon[0] = ucon_norm(ucon, gcov); + + // now get three velocity + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real Wd = lapsed * ucon[0]; + Real W = lapse * ucon[0]; - // finally compute three-velocity for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + Wd * beta[d] / lapsed; + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } + printf("POSTSHOCK r = %g u^0 (should match vr1) = %g\n", r, v(ivlo, k, j, i)); + + // preshock - 0 } else { - // postshock - 1 + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; - Real vr0 = -1. * std::sqrt(W0 * W0 - 1.) / W0; + Real vr0 = std::sqrt(lapse0 * lapse0 - 1.); Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - Real alphasq = 1. - (2. / r); - Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); - Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - Real rho1 = rho0 * W0 * (vr0 / vr1); + printf("PRESHOCK: passing values of: r, rho0, vr0 = %g %g %g\n", r, rho0, vr0); + Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, + eos_lambda[0]); - const Real epsND = 0.003; // 2.7e16 ergs / g for M = 1.3Mpns + // printf("PRESHOCK: r = %e T (K) = %e\n", r, T/pc.kb); - // printf("postshock r= 0, "discr < 0"); + return (-BB - std::sqrt(discr)) / (2. * AA); +} + } // namespace standing_accretion_shock From b60b96d4c97997b941c849aa6d113fe5e1dc1130 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 8 Jun 2023 14:24:39 -0600 Subject: [PATCH 30/51] cleanup, remove prints, this runs and completes setup then fails to bracket root when executing driver --- src/pgen/standing_accretion_shock.cpp | 32 +-------------------------- 1 file changed, 1 insertion(+), 31 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 81f52ff46..c407a0e94 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -39,8 +39,6 @@ class MachResidual { Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity - // printf("Mach Residual Func produced: u, w, cs, Mach, vr0 = %g %g %g %g %g\n", u, w, - // cs, mach, vr0_); return mach - target_mach_; } @@ -50,7 +48,6 @@ class MachResidual { Real lambda_[2]; }; -// protehypes---- KOKKOS_FUNCTION Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, const Real target_mach, const Real Tmin, const Real Tmax, @@ -58,8 +55,6 @@ Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); -//-------------- - void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE( @@ -131,10 +126,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real gamma = 4. / 3.; const Real epsND = 0.003; // 2.7e16 ergs / g for M = 1.3Mpns - // printf("0 (preshock region) x1, r, rho0, vr0, lapse0, W0 = %g %g %g %g %g - // %g\n",x1, r, rho0, vr0, lapse0, W0); - - // set Ye everywhere Real eos_lambda[2]; if (iye > 0) { v(iye, k, j, i) = 0.5; @@ -154,10 +145,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real vr1 = -1. * (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); - printf("inside shock (postshock region) r, rho1, vr1 = %g %g %g\n", r, rho1, - vr1); - // printf("inside shock (postshock region) r, alphasq, psi, rho0, vr0, rho1, - // vr1 = %g %g %g %g %g %g %g\n",r, alphasq, psi, rho0, vr0, rho1, vr1); v(irho, k, j, i) = rho1; v(ieng, k, j, i) = (W0 - 1. + epsND * (gamma - 1.)) / (gamma); v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( @@ -168,15 +155,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - // construct contravariant 4-velocity Real ucon[] = {0.0, vr1, 0.0, 0.0}; - - // initialize spacetime metric Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); ucon[0] = ucon_norm(ucon, gcov); - // now get three velocity const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); @@ -186,22 +169,17 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } - printf("POSTSHOCK r = %g u^0 (should match vr1) = %g\n", r, v(ivlo, k, j, i)); - // preshock - 0 } else { Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); Real W0 = 1. / lapse0; - Real vr0 = std::sqrt(lapse0 * lapse0 - 1.); + Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - printf("PRESHOCK: passing values of: r, rho0, vr0 = %g %g %g\n", r, rho0, vr0); Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, eos_lambda[0]); - // printf("PRESHOCK: r = %e T (K) = %e\n", r, T/pc.kb); - v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T; v(ieng, k, j, i) = @@ -212,15 +190,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - // construct contravariant 4-velocity Real ucon[] = {0.0, vr0, 0.0, 0.0}; - - // initialize spacetime metric Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); ucon[0] = ucon_norm(ucon, gcov); - // now get three velocity const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); Real beta[3]; geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); @@ -229,8 +203,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { for (int d = 0; d < 3; d++) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } - - printf("preshock r = %g u^0 (should match vr0) = %g\n", r, v(ivlo, k, j, i)); } }); @@ -241,8 +213,6 @@ KOKKOS_FUNCTION Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, const Real Tmin, const Real Tmax, const Real vr0, const Real Ye) { - // printf("passing values of: rho, vr0, Mach_target, Ye = %g %g %g %g\n", rho, vr0, - // target_mach, Ye); MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; Real Troot = From 626de315f58a781c6537b53faa6aebd92af22de6 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 14 Jun 2023 10:29:32 -0600 Subject: [PATCH 31/51] testing --- src/pgen/standing_accretion_shock.cpp | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index c407a0e94..038a6b2a9 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -16,6 +16,7 @@ #include "phoebus_utils/root_find.hpp" #include "phoebus_utils/unit_conversions.hpp" #include "utils/error_checking.hpp" +#include namespace standing_accretion_shock { @@ -135,14 +136,14 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // postshock - 1 if (r < rShock) { - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); - Real W0 = 1. / lapse0; - Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); + const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real alphasq = 1. - (2. / r); Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); - Real vr1 = -1. * (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); v(irho, k, j, i) = rho1; @@ -155,6 +156,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g " + "%g %g %g\n", + r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + v(iprs, k, j, i), v(igm1, k, j, i)); + Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); @@ -172,10 +178,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // preshock - 0 } else { - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - Real W0 = 1. / lapse0; + Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); - Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); + Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, eos_lambda[0]); @@ -190,6 +196,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g " + "%g\n", + r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + v(iprs, k, j, i), v(igm1, k, j, i)); + Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); From e47dd1798df9c34602d08970f66696f22897fc73 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 14 Jun 2023 14:35:16 -0600 Subject: [PATCH 32/51] fix rs conversion, try solving T for a target mach to close postshock vars --- inputs/standing_accretion_shock.pin | 17 +++---- src/pgen/standing_accretion_shock.cpp | 71 +++++++++++---------------- 2 files changed, 38 insertions(+), 50 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 5913479ef..52117474b 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -12,7 +12,7 @@ problem = standing_accretion_shock -ix1_bc = outflow +ix1_bc = reflect ox1_bc = outflow ix2_bc = reflect ox2_bc = reflect @@ -29,18 +29,19 @@ variables = p.density, & c.energy, & pressure, & cs, & + temperature,& p.ye, & g.c.coord, & g.n.coord file_type = hdf5 # Tabular data dump -dt = 1.e-3 # time increment between outputs +dt = 100 # time increment between outputs nlim = -1 # cycle limit -tlim = 1.0 # time limit +tlim = 5000 # time limit integrator = rk2 # time integration algorithm -ncycle_out = 1 # interval for stdout summary info +ncycle_out = 100 # interval for stdout summary info dt_init_fact = 1.e-6 @@ -49,8 +50,8 @@ nghost = 4 #numlevel = 3 nx1 = 256 # Number of zones in X1-direction -x1min = 3.44199 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 -x1max = 5.5 # maximum value of X1 | rmax | x1max => 472.5 (km) +x1min = 31.2942 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 +x1max = 246.087 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag @@ -103,9 +104,6 @@ mhd = false a = 0 - -derefine_poles = false - scale_free = false geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km @@ -115,3 +113,4 @@ fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km target_mach = 100 # target Mach number for upstream flow in preshock region +target_mach_ps = 1 # target Mach number for postshock region diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 038a6b2a9..7866b9c81 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -88,8 +88,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - const Real target_mach = - pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); + const Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); + const Real target_mach_ps = pin->GetOrAddReal("standing_accretion_shock", "target_mach_ps", 1); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -108,9 +108,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - const Real MPNS = 1.3 * solar_mass * unit_conv.GetMassCGSToCode(); - const Real rs = - (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + const Real MPNS = 1.3 * solar_mass; // grams + const Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); printf("rShock (code units) = %g \n", rShock); @@ -139,27 +138,25 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + const Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - Real alphasq = 1. - (2. / r); - Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); - Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + //Real alphasq = 1. - (rs / r); + const Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); + Real vr1 = -1. * std::sqrt(lapse1 * lapse1 - 1.); + + //Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); + //Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); - v(irho, k, j, i) = rho1; - v(ieng, k, j, i) = (W0 - 1. + epsND * (gamma - 1.)) / (gamma); - v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy( - rho1, v(ieng, k, j, i) / rho1, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); - - printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g " - "%g %g %g\n", - r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - v(iprs, k, j, i), v(igm1, k, j, i)); + Real T1 = temperature_from_rho_mach(eos, rho1, target_mach_ps, Tmin, Tmax, vr1, eos_lambda[0]); + + v(irho, k, j, i) = rho1; + v(itmp, k, j, i) = T1; + v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + + printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -175,31 +172,23 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } - // preshock - 0 + // preshock - 0 } else { - Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); - Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); + const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + const Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, - eos_lambda[0]); + Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, eos_lambda[0]); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T; - v(ieng, k, j, i) = - rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); - - printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g " - "%g\n", - r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - v(iprs, k, j, i), v(igm1, k, j, i)); + v(ieng, k, j, i) = rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); + + printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; From b9a0d26b8baf1ab2d14a28936518366c3a157906 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 14 Jun 2023 15:10:25 -0600 Subject: [PATCH 33/51] actually try this --- src/pgen/standing_accretion_shock.cpp | 66 +++++++++++++++++---------- 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 7866b9c81..9deaf3b84 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -88,8 +88,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - const Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); - const Real target_mach_ps = pin->GetOrAddReal("standing_accretion_shock", "target_mach_ps", 1); + const Real target_mach = + pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); + const Real target_mach_ps = + pin->GetOrAddReal("standing_accretion_shock", "target_mach_ps", 1); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -109,7 +111,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); const Real MPNS = 1.3 * solar_mass; // grams - const Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + const Real rs = + (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); printf("rShock (code units) = %g \n", rShock); @@ -139,24 +142,29 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); const Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - - //Real alphasq = 1. - (rs / r); + const Real alphasq = 1. - (rs / r); const Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); - Real vr1 = -1. * std::sqrt(lapse1 * lapse1 - 1.); + const Real vr1 = -1. * std::sqrt(lapse1 * lapse1 - 1.); - //Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); - //Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; Real rho1 = rho0 * W0 * (vr0 / vr1); - - Real T1 = temperature_from_rho_mach(eos, rho1, target_mach_ps, Tmin, Tmax, vr1, eos_lambda[0]); - - v(irho, k, j, i) = rho1; - v(itmp, k, j, i) = T1; - v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature(rho1, T1, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - - printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), v(iprs, k, j, i), v(igm1, k, j, i)); + Real h1 = (rho0 * vr0 * W0 * W0) / (rho1 * vr1); + Real p1 = (rho0 * vr0 * vr0 * W0 * W0 - rho0 * vr0 * vr1 * W0 * W0) / alphasq; + Real eps1 = h1 - 1. - p1 / rho1; + + v(irho, k, j, i) = rho1; + v(itmp, k, j, i) = + eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); + v(ieng, k, j, i) = rho1 * eps1; + v(iprs, k, j, i) = eos.PressureFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); + + printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g " + "%g %g %g\n", + r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -172,7 +180,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } - // preshock - 0 + // preshock - 0 } else { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); @@ -180,15 +188,23 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); const Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, eos_lambda[0]); + Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, + eos_lambda[0]); v(irho, k, j, i) = rho0; v(itmp, k, j, i) = T; - v(ieng, k, j, i) = rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); - v(iprs, k, j, i) = eos.PressureFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - - printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); + v(ieng, k, j, i) = + rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); + v(iprs, k, j, i) = eos.PressureFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); + + printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g " + "%g\n", + r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; From a5956e20eb2c9c8888653da3cd5cc9df5375a9fc Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Fri, 16 Jun 2023 14:30:36 -0600 Subject: [PATCH 34/51] try this, low temp cell somewhere causes failure in root find --- src/pgen/standing_accretion_shock.cpp | 33 ++++++++++++++------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 9deaf3b84..0a9d1491c 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -127,7 +127,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real r = std::abs(x1); const Real gamma = 4. / 3.; - const Real epsND = 0.003; // 2.7e16 ergs / g for M = 1.3Mpns Real eos_lambda[2]; if (iye > 0) { @@ -141,30 +140,29 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + const Real alphasq = 1. - (rs / r); - const Real lapse1 = geom.Lapse(CellLocation::Cent, k, j, r); - const Real vr1 = -1. * std::sqrt(lapse1 * lapse1 - 1.); + const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); + const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + const Real rho1 = rho0 * W0 * (vr0 / vr1); - Real rho1 = rho0 * W0 * (vr0 / vr1); Real h1 = (rho0 * vr0 * W0 * W0) / (rho1 * vr1); - Real p1 = (rho0 * vr0 * vr0 * W0 * W0 - rho0 * vr0 * vr1 * W0 * W0) / alphasq; + Real p1 = (-h1 * rho1 * vr1 * vr1 + rho0 * vr0 * vr0 * W0 * W0) / alphasq; Real eps1 = h1 - 1. - p1 / rho1; v(irho, k, j, i) = rho1; v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); v(ieng, k, j, i) = rho1 * eps1; - v(iprs, k, j, i) = eos.PressureFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); + v(iprs, k, j, i) = p1; v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - printf("POSTSHOCK r, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g " - "%g %g %g\n", - r, W0, vr1, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - v(iprs, k, j, i), v(igm1, k, j, i)); + // printf("POSTSHOCK r, alphasq, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g + // %g %g %g %g %g %g %g\n",r, alphasq,W0, vr1, v(irho, k, j, i), v(ieng, k, j, + // i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -201,10 +199,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g " - "%g\n", - r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - v(iprs, k, j, i), v(igm1, k, j, i)); + // printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g + // %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, + // i),v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; @@ -220,6 +217,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } } + + printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r, + v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + v(iprs, k, j, i), v(igm1, k, j, i)); }); fluid::PrimitiveToConserved(rc); From c04b29cc0faf000bc5d7791abd42d4ff441ad5ad Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 20 Jun 2023 09:27:39 -0600 Subject: [PATCH 35/51] debugging --- src/pgen/standing_accretion_shock.cpp | 49 +++++++++++++++++---------- 1 file changed, 32 insertions(+), 17 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 0a9d1491c..a40912459 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -11,6 +11,10 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. + +// Parthenon +#include + #include "geometry/boyer_lindquist.hpp" #include "pgen/pgen.hpp" #include "phoebus_utils/root_find.hpp" @@ -40,7 +44,8 @@ class MachResidual { Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity - return mach - target_mach_; + //printf("rho, vr0, target_mach, P, eps, cs, mach, T = %g %g %g %g %g %g %g %g\n", rho_, vr0_, target_mach_, P, eps, cs, mach, T); + return target_mach_ - mach; } private: @@ -88,10 +93,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - const Real target_mach = - pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); - const Real target_mach_ps = - pin->GetOrAddReal("standing_accretion_shock", "target_mach_ps", 1); + Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -101,18 +103,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); - const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - auto &unit_conv = - pmb->packages.Get("phoebus")->Param("unit_conv"); + auto &unit_conv = pmb->packages.Get("phoebus")->Param("unit_conv"); Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - const Real MPNS = 1.3 * solar_mass; // grams - const Real rs = - (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + Real MPNS = 1.3 * solar_mass; + Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); printf("rShock (code units) = %g \n", rShock); @@ -151,6 +150,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real p1 = (-h1 * rho1 * vr1 * vr1 + rho0 * vr0 * vr0 * W0 * W0) / alphasq; Real eps1 = h1 - 1. - p1 / rho1; + if (std::isnan(vr0) || std::isnan(rho0) || std::isnan(vr1) || std::isnan(rho1)) { + printf("std::isnan(vr0) || std::isnan(rho0) || std::isnan(vr1) || std::isnan(rho1) @ r = %i %i %i %i %g \n", std::isnan(vr0), std::isnan(rho0), std::isnan(vr1), std::isnan(rho1), r); + + } + v(irho, k, j, i) = rho1; v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); @@ -181,10 +185,17 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // preshock - 0 } else { - const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, r); + const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + const Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + //printf("r, lapse0, vr0, W0, rho0 = %g %g %g %g %g \n", r, lapse0, vr0, W0, rho0); + + if (std::isnan(vr0) || std::isnan(rho0) ) { + printf("std::isnan(vr0) || std::isnan(rho0) @ r = %i %i %g \n", std::isnan(vr0), std::isnan(rho0), r); + + } + Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, eos_lambda[0]); @@ -218,9 +229,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } } - printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r, - v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - v(iprs, k, j, i), v(igm1, k, j, i)); + //int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); + //PARTHENON_REQUIRE(num_nans == 0, "NaNs founds"); + //printf("num_nans, r = %i %g \n", num_nans, r); + + //printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r, + // v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), + // v(iprs, k, j, i), v(igm1, k, j, i)); }); fluid::PrimitiveToConserved(rc); @@ -233,7 +248,7 @@ Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; Real Troot = - root.regula_falsi(res, Tmin, Tmax, 1.e-10 * target_mach, std::max(Tmin, 1.e-10)); + root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, std::max(Tmin, 1.e-3)); return Troot; } From 8b1d6c28a718f978751999b5cf176506e8b42fc3 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 21 Jun 2023 10:23:00 -0600 Subject: [PATCH 36/51] assumes ideal gas for vel --- src/pgen/standing_accretion_shock.cpp | 71 +++++++++++++-------------- 1 file changed, 33 insertions(+), 38 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index a40912459..e03bd9e70 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -45,7 +45,7 @@ class MachResidual { Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity //printf("rho, vr0, target_mach, P, eps, cs, mach, T = %g %g %g %g %g %g %g %g\n", rho_, vr0_, target_mach_, P, eps, cs, mach, T); - return target_mach_ - mach; + return mach - target_mach_; } private: @@ -137,36 +137,36 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (r < rShock) { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + const Real vr0 = -1.*std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - - const Real alphasq = 1. - (rs / r); + const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); + const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real T0 = temperature_from_rho_mach(eos, rho_0_Shock, target_mach, Tmin, Tmax, vr0,eos_lambda[0]); + const Real p0 = eos.PressureFromDensityTemperature(rho_0_Shock, T0, eos_lambda); + const Real eps0 = eos.InternalEnergyFromDensityTemperature(rho_0_Shock, T0, eos_lambda); + const Real h0 = 1. + eps0 + p0/rho_0_Shock; + + const Real alphasq = 1. - (rs / r); const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - const Real rho1 = rho0 * W0 * (vr0 / vr1); - - Real h1 = (rho0 * vr0 * W0 * W0) / (rho1 * vr1); - Real p1 = (-h1 * rho1 * vr1 * vr1 + rho0 * vr0 * vr0 * W0 * W0) / alphasq; - Real eps1 = h1 - 1. - p1 / rho1; + const Real W1 = 1. / sqrt(1. - std::pow(vr1, 2)); + const Real rho1 = (rho_0 * W0 * vr0) / (W1 * vr1); + const Real h1 = (h0 * vr0 * W0 * W0 * rho_0) / (W1 * W1 * rho1 * vr1); + const Real p1 = (alphasq * p0 + h0 * vr0 * vr0 * W0 * W0 * rho_0 - h1 * vr1 * vr1 * W1 * W1 * rho1) / alphasq; + const Real eps1 = h1 - 1. - p1 / rho1; - if (std::isnan(vr0) || std::isnan(rho0) || std::isnan(vr1) || std::isnan(rho1)) { - printf("std::isnan(vr0) || std::isnan(rho0) || std::isnan(vr1) || std::isnan(rho1) @ r = %i %i %i %i %g \n", std::isnan(vr0), std::isnan(rho0), std::isnan(vr1), std::isnan(rho1), r); + if (std::isnan(vr0) || std::isnan(rho_0) || std::isnan(vr1) || std::isnan(rho1)) { + printf("std::isnan(vr0) || std::isnan(rho_0) || std::isnan(vr1) || std::isnan(rho1) @ r = %i %i %i %i %g \n", std::isnan(vr0), std::isnan(rho_0), std::isnan(vr1), std::isnan(rho1), r); } v(irho, k, j, i) = rho1; - v(itmp, k, j, i) = - eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); + v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); v(ieng, k, j, i) = rho1 * eps1; v(iprs, k, j, i) = p1; - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - // printf("POSTSHOCK r, alphasq, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g - // %g %g %g %g %g %g %g\n",r, alphasq,W0, vr1, v(irho, k, j, i), v(ieng, k, j, - // i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); + //printf("POSTSHOCK r, alphasq, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g %g %g %g %g\n",r, alphasq,W0, vr1, v(irho, k, j, i), v(ieng, k, j,i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -186,33 +186,30 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } else { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); + const Real vr0 = -1.*std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - //printf("r, lapse0, vr0, W0, rho0 = %g %g %g %g %g \n", r, lapse0, vr0, W0, rho0); + const Real rho_0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + //printf("r, lapse0, vr0, W0, rho_0 = %g %g %g %g %g \n", r, lapse0, vr0, W0, rho_0); - if (std::isnan(vr0) || std::isnan(rho0) ) { - printf("std::isnan(vr0) || std::isnan(rho0) @ r = %i %i %g \n", std::isnan(vr0), std::isnan(rho0), r); + if (std::isnan(vr0) || std::isnan(rho_0) ) { + printf("std::isnan(vr0) || std::isnan(rho_0) @ r = %i %i %g \n", std::isnan(vr0), std::isnan(rho_0), r); } - Real T = temperature_from_rho_mach(eos, rho0, target_mach, Tmin, Tmax, vr0, - eos_lambda[0]); + Real T = temperature_from_rho_mach(eos, rho_0, target_mach, Tmin, Tmax, vr0,eos_lambda[0]); - v(irho, k, j, i) = rho0; + v(irho, k, j, i) = rho_0; v(itmp, k, j, i) = T; v(ieng, k, j, i) = - rho0 * eos.InternalEnergyFromDensityTemperature(rho0, T, eos_lambda); + rho_0 * eos.InternalEnergyFromDensityTemperature(rho_0, T, eos_lambda); v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - // printf("PRESHOCK r, W0, vr0, rho0, eint, T0, P0, gma1 = %g %g %g %g %g %g - // %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, - // i),v(iprs, k, j, i), v(igm1, k, j, i)); + //printf("PRESHOCK r, W0, vr0, rho_0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; @@ -229,13 +226,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } } - //int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); - //PARTHENON_REQUIRE(num_nans == 0, "NaNs founds"); + int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); + PARTHENON_REQUIRE(num_nans == 0, "NaNs founds"); //printf("num_nans, r = %i %g \n", num_nans, r); - //printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r, - // v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i), - // v(iprs, k, j, i), v(igm1, k, j, i)); + //printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r,v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); }); fluid::PrimitiveToConserved(rc); @@ -248,7 +243,7 @@ Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; Real Troot = - root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, std::max(Tmin, 1.e-3)); + root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, Tmin-1e-10); return Troot; } From 9f021a31a18b83609aa4d1bcfbd345887407822e Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 22 Jun 2023 16:02:51 -0600 Subject: [PATCH 37/51] doesnt work --- inputs/standing_accretion_shock.pin | 18 +++- src/pgen/standing_accretion_shock.cpp | 129 ++++++++++++++++---------- 2 files changed, 95 insertions(+), 52 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 52117474b..68bd7c7bb 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -49,7 +49,7 @@ nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 256 # Number of zones in X1-direction +nx1 = 128 # Number of zones in X1-direction x1min = 31.2942 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 x1max = 246.087 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag @@ -57,7 +57,7 @@ ox1_bc = user # Outer-X1 boundary condition flag nx2 = 1 # Number of zones in X2-direction x2min = 0 # minimum value of X2 -x2max = 1 #3.14159265359 #1 # maximum value of X2. Pi +x2max = 3.14159265359 #3.14159265359 #1 # maximum value of X2. Pi for BL ix2_bc = user # Inner-X2 boundary condition flag ox2_bc = user # Outer-X2 boundary condition flag @@ -101,6 +101,17 @@ c2p_tol = 1.e-8 Ye = true mhd = false + +enable_floors = true +enable_ceilings = true +enable_flux_fixup = true +enable_c2p_fixup = true +floor_type = ConstantRhoSie +rho0_floor = 1.e-2 +sie0_floor = 1.e4 +rho0_ceiling = 2. +sie0_ceiling = 1.e7 + a = 0 @@ -113,4 +124,5 @@ fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km target_mach = 100 # target Mach number for upstream flow in preshock region -target_mach_ps = 1 # target Mach number for postshock region +rhomin = 1.e-2 +rhomax = 3. diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index e03bd9e70..80485510a 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -11,7 +11,6 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. - // Parthenon #include @@ -44,7 +43,6 @@ class MachResidual { Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity - //printf("rho, vr0, target_mach, P, eps, cs, mach, T = %g %g %g %g %g %g %g %g\n", rho_, vr0_, target_mach_, P, eps, cs, mach, T); return mach - target_mach_; } @@ -54,6 +52,29 @@ class MachResidual { Real lambda_[2]; }; +class EnergyResidual { + public: + KOKKOS_FUNCTION + EnergyResidual(const Real alphasq, const Real p0, const Real h0, const Real v0, + const Real W0, const Real rho0, const Real W1) + : alphasq_(alphasq), p0_(p0), h0_(h0), v0_(v0), W0_(W0), rho0_(rho0), W1_(W1) {} + KOKKOS_INLINE_FUNCTION + Real operator()(const Real rho1) { + Real p1 = (-h0_ * std::pow(v0_, 2) * std::pow(W0_, 3) * std::pow(rho0_, 2) + + alphasq_ * p0_ * W1_ * rho1 + + h0_ * std::pow(v0_, 2) * std::pow(W0_, 2) * W1_ * rho0_ * rho1) / + (alphasq_ * W1_ * rho1); + Real v1 = (v0_ * W0_ * rho0_) / (W1_ * rho1); + Real h1 = (h0_ * v0_ * std::pow(W0_, 2) * rho0_) / (v1 * std::pow(W1_, 2) * rho1); + Real lhs = rho1 * h1 * std::pow(W1_, 2) * std::pow(v1, 2) + p1 * alphasq_; + Real rhs = rho0_ * h0_ * std::pow(W0_, 2) * std::pow(v0_, 2) + p0_ * alphasq_; + return lhs - rhs; + } + + private: + Real alphasq_, p0_, h0_, v0_, W0_, rho0_, W1_; +}; + KOKKOS_FUNCTION Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, const Real target_mach, const Real Tmin, const Real Tmax, @@ -61,6 +82,11 @@ Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); +KOKKOS_FUNCTION +Real density_from_jump_conditions(const Real alphasq, const Real p0, const Real h0, + const Real v0, const Real W0, const Real rho0, + const Real W1, const Real Pmin, const Real Pmax); + void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE( @@ -94,6 +120,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); + Real rhomin = pin->GetOrAddReal("standing_accretion_shock", "rhomin", 1e-2); + Real rhomax = pin->GetOrAddReal("standing_accretion_shock", "rhomax", 3); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); @@ -105,13 +133,14 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto Tmax = pmb->packages.Get("eos")->Param("T_max"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - - auto &unit_conv = pmb->packages.Get("phoebus")->Param("unit_conv"); + auto &unit_conv = + pmb->packages.Get("phoebus")->Param("unit_conv"); Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); Real MPNS = 1.3 * solar_mass; - Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + Real rs = + (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); printf("rShock (code units) = %g \n", rShock); @@ -125,8 +154,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x3 = coords.Xc<3>(k, j, i); Real r = std::abs(x1); - const Real gamma = 4. / 3.; - Real eos_lambda[2]; if (iye > 0) { v(iye, k, j, i) = 0.5; @@ -137,36 +164,40 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (r < rShock) { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1.*std::sqrt(lapse0 * lapse0 - 1.); + const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - const Real T0 = temperature_from_rho_mach(eos, rho_0_Shock, target_mach, Tmin, Tmax, vr0,eos_lambda[0]); - const Real p0 = eos.PressureFromDensityTemperature(rho_0_Shock, T0, eos_lambda); - const Real eps0 = eos.InternalEnergyFromDensityTemperature(rho_0_Shock, T0, eos_lambda); - const Real h0 = 1. + eps0 + p0/rho_0_Shock; - - const Real alphasq = 1. - (rs / r); - const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1.) / W0); - const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - const Real W1 = 1. / sqrt(1. - std::pow(vr1, 2)); - const Real rho1 = (rho_0 * W0 * vr0) / (W1 * vr1); - const Real h1 = (h0 * vr0 * W0 * W0 * rho_0) / (W1 * W1 * rho1 * vr1); - const Real p1 = (alphasq * p0 + h0 * vr0 * vr0 * W0 * W0 * rho_0 - h1 * vr1 * vr1 * W1 * W1 * rho1) / alphasq; + const Real rho_0_Shock = + Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real T0 = temperature_from_rho_mach(eos, rho_0_Shock, target_mach, Tmin, + Tmax, vr0, eos_lambda[0]); + const Real p0 = eos.PressureFromDensityTemperature(rho_0_Shock, T0, eos_lambda); + const Real eps0 = + eos.InternalEnergyFromDensityTemperature(rho_0_Shock, T0, eos_lambda); + const Real h0 = 1. + eps0 + p0 / rho_0_Shock; + const Real alphasq = 1. - (rs / r); + const Real W1 = 1.; + + Real rho1 = density_from_jump_conditions(alphasq, p0, h0, vr0, W0, rho_0_Shock, + W1, rhomin, rhomax); + const Real p1 = + (-h0 * std::pow(vr0, 2) * std::pow(W0, 3) * std::pow(rho_0_Shock, 2) + + alphasq * p0 * W1 * rho1 + + h0 * std::pow(vr0, 2) * std::pow(W0, 2) * W1 * rho_0_Shock * rho1) / + (alphasq * W1 * rho1); + const Real vr1 = (vr0 * W0 * rho_0_Shock) / (W1 * rho1); + const Real h1 = + (h0 * vr0 * std::pow(W0, 2) * rho_0_Shock) / (vr1 * std::pow(W1, 2) * rho1); const Real eps1 = h1 - 1. - p1 / rho1; - if (std::isnan(vr0) || std::isnan(rho_0) || std::isnan(vr1) || std::isnan(rho1)) { - printf("std::isnan(vr0) || std::isnan(rho_0) || std::isnan(vr1) || std::isnan(rho1) @ r = %i %i %i %i %g \n", std::isnan(vr0), std::isnan(rho_0), std::isnan(vr1), std::isnan(rho1), r); - - } - v(irho, k, j, i) = rho1; - v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); + v(itmp, k, j, i) = + eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); v(ieng, k, j, i) = rho1 * eps1; v(iprs, k, j, i) = p1; - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - - //printf("POSTSHOCK r, alphasq, W0,vr1, rho1, eint1, T1(K), P1, gma1 = %g %g %g %g %g %g %g %g %g\n",r, alphasq,W0, vr1, v(irho, k, j, i), v(ieng, k, j,i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( + v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / + v(iprs, k, j, i); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -186,18 +217,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } else { const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1.*std::sqrt(lapse0 * lapse0 - 1.); + const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho_0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - //printf("r, lapse0, vr0, W0, rho_0 = %g %g %g %g %g \n", r, lapse0, vr0, W0, rho_0); - - if (std::isnan(vr0) || std::isnan(rho_0) ) { - printf("std::isnan(vr0) || std::isnan(rho_0) @ r = %i %i %g \n", std::isnan(vr0), std::isnan(rho_0), r); - - } - - - Real T = temperature_from_rho_mach(eos, rho_0, target_mach, Tmin, Tmax, vr0,eos_lambda[0]); + const Real rho_0 = + Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + Real T = temperature_from_rho_mach(eos, rho_0, target_mach, Tmin, Tmax, vr0, + eos_lambda[0]); v(irho, k, j, i) = rho_0; v(itmp, k, j, i) = T; @@ -209,8 +234,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - //printf("PRESHOCK r, W0, vr0, rho_0, eint, T0, P0, gma1 = %g %g %g %g %g %g %g %g\n",r, W0, vr0, v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); - Real ucon[] = {0.0, vr0, 0.0, 0.0}; Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); @@ -226,11 +249,10 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { } } - int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); + int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + + std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + + std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); PARTHENON_REQUIRE(num_nans == 0, "NaNs founds"); - //printf("num_nans, r = %i %g \n", num_nans, r); - - //printf("r, vr, rho, eint, T, P, gma1 = %g %g %g %g %g %g %g\n", r,v(ivlo, k, j, i), v(irho, k, j, i), v(ieng, k, j, i), v(itmp, k, j, i),v(iprs, k, j, i), v(igm1, k, j, i)); }); fluid::PrimitiveToConserved(rc); @@ -242,8 +264,7 @@ Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target const Real Ye) { MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; - Real Troot = - root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, Tmin-1e-10); + Real Troot = root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, Tmin - 1e-10); return Troot; } @@ -261,4 +282,14 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]) { return (-BB - std::sqrt(discr)) / (2. * AA); } +KOKKOS_FUNCTION +Real density_from_jump_conditions(const Real alphasq, const Real p0, const Real h0, + const Real v0, const Real W0, const Real rho0, + const Real W1, const Real rhomin, const Real rhomax) { + EnergyResidual res(alphasq, p0, h0, v0, W0, rho0, W1); + root_find::RootFind root; + Real rho1root = root.regula_falsi(res, rhomin, rhomax, 1.e-10, rhomin * 5); + return rho1root; +} + } // namespace standing_accretion_shock From 30dfc1a6bb3246144888ec052e674bed420be3fc Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 27 Jun 2023 15:58:02 -0600 Subject: [PATCH 38/51] ideal gas, not working --- inputs/standing_accretion_shock.pin | 42 ++++--- src/pgen/standing_accretion_shock.cpp | 154 ++++++++------------------ 2 files changed, 68 insertions(+), 128 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 68bd7c7bb..b91d9e160 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -35,21 +35,21 @@ variables = p.density, & g.n.coord file_type = hdf5 # Tabular data dump -dt = 100 # time increment between outputs +dt = 5.e-3 # time increment between outputs -nlim = -1 # cycle limit +nlim = 100 # cycle limit tlim = 5000 # time limit -integrator = rk2 # time integration algorithm -ncycle_out = 100 # interval for stdout summary info -dt_init_fact = 1.e-6 +integrator = rk4 # time integration algorithm +ncycle_out = 1 # interval for stdout summary info +dt_init_fact = 1.e-9 nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 128 # Number of zones in X1-direction +nx1 = 256 # Number of zones in X1-direction x1min = 31.2942 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 x1max = 246.087 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag @@ -80,9 +80,9 @@ method = derivative_order_1 max_level = 3 -type = StellarCollapse -filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 -use_sp5 = false +type = IdealGas +Gamma = 1.3333333333333333 +Cv = 3.0 hydro = true @@ -92,25 +92,23 @@ rad = false xorder = 2 -cfl = 0.8 +cfl = 0.5 riemann = llf -recon = linear +recon = weno5 c2p_method = robust -c2p_max_iter = 100 +c2p_max_iter = 1000 c2p_tol = 1.e-8 Ye = true mhd = false -enable_floors = true -enable_ceilings = true -enable_flux_fixup = true -enable_c2p_fixup = true +enable_floors = false +enable_ceilings = false +enable_flux_fixup = false +enable_c2p_fixup = false floor_type = ConstantRhoSie -rho0_floor = 1.e-2 -sie0_floor = 1.e4 -rho0_ceiling = 2. -sie0_ceiling = 1.e7 +rho0_floor = 1.e-3 +sie0_floor = 1.e-8 a = 0 @@ -123,6 +121,4 @@ fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km -target_mach = 100 # target Mach number for upstream flow in preshock region -rhomin = 1.e-2 -rhomax = 3. +target_mach = -100 # target Mach number for upstream flow in preshock region diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 80485510a..0234d35fc 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -35,14 +35,10 @@ class MachResidual { lambda_[0] = Ye; } KOKKOS_INLINE_FUNCTION - Real operator()(const Real T) { - Real P = eos_.PressureFromDensityTemperature(rho_, T, lambda_); - Real eps = eos_.InternalEnergyFromDensityTemperature(rho_, T, lambda_); - Real bmod = eos_.BulkModulusFromDensityTemperature(rho_, T, lambda_); - Real u = rho_ * eps; // convert eps / V to specific internal energy - Real w = rho_ + P + u; // h = 1 + eps + P/rho | w = rho * h == rho + u + P - Real cs = std::sqrt(bmod / w); // cs^2 = bmod / w - Real mach = std::abs(vr0_) / cs; // radial component of preshock velocity + Real operator()(const Real eps) { + const Real gamma = 4. / 3.; + Real cs = std::sqrt(gamma * (gamma - 1.) * eps); + Real mach = vr0_ / cs; return mach - target_mach_; } @@ -52,41 +48,13 @@ class MachResidual { Real lambda_[2]; }; -class EnergyResidual { - public: - KOKKOS_FUNCTION - EnergyResidual(const Real alphasq, const Real p0, const Real h0, const Real v0, - const Real W0, const Real rho0, const Real W1) - : alphasq_(alphasq), p0_(p0), h0_(h0), v0_(v0), W0_(W0), rho0_(rho0), W1_(W1) {} - KOKKOS_INLINE_FUNCTION - Real operator()(const Real rho1) { - Real p1 = (-h0_ * std::pow(v0_, 2) * std::pow(W0_, 3) * std::pow(rho0_, 2) + - alphasq_ * p0_ * W1_ * rho1 + - h0_ * std::pow(v0_, 2) * std::pow(W0_, 2) * W1_ * rho0_ * rho1) / - (alphasq_ * W1_ * rho1); - Real v1 = (v0_ * W0_ * rho0_) / (W1_ * rho1); - Real h1 = (h0_ * v0_ * std::pow(W0_, 2) * rho0_) / (v1 * std::pow(W1_, 2) * rho1); - Real lhs = rho1 * h1 * std::pow(W1_, 2) * std::pow(v1, 2) + p1 * alphasq_; - Real rhs = rho0_ * h0_ * std::pow(W0_, 2) * std::pow(v0_, 2) + p0_ * alphasq_; - return lhs - rhs; - } - - private: - Real alphasq_, p0_, h0_, v0_, W0_, rho0_, W1_; -}; - KOKKOS_FUNCTION -Real temperature_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, - const Real target_mach, const Real Tmin, const Real Tmax, - const Real vr0, const Real Ye); +Real eps_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, + const Real target_mach, const Real epsmin, const Real epsmax, + const Real vr0, const Real Ye); KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); -KOKKOS_FUNCTION -Real density_from_jump_conditions(const Real alphasq, const Real p0, const Real h0, - const Real v0, const Real W0, const Real rho0, - const Real W1, const Real Pmin, const Real Pmax); - void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { PARTHENON_REQUIRE( @@ -114,27 +82,26 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const std::string eos_type = pin->GetString("eos", "type"); PARTHENON_REQUIRE_THROWS( - eos_type == "StellarCollapse", - "Standing Accretion Shock setup only works with StellarCollapse EOS"); + eos_type == "IdealGas", + "Standing Accretion Shock setup only works with Ideal Gas EOS"); Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", 100); - Real rhomin = pin->GetOrAddReal("standing_accretion_shock", "rhomin", 1e-2); - Real rhomax = pin->GetOrAddReal("standing_accretion_shock", "rhomax", 3); + Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", -100); IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); auto &coords = pmb->coords; + auto &unit_conv = + pmb->packages.Get("phoebus")->Param("unit_conv"); auto eos = pmb->packages.Get("eos")->Param("d.EOS"); - auto Tmin = pmb->packages.Get("eos")->Param("T_min"); - auto Tmax = pmb->packages.Get("eos")->Param("T_max"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - auto &unit_conv = - pmb->packages.Get("phoebus")->Param("unit_conv"); + + const Real epsmin = 1.e-9; + const Real epsmax = 1.e6; Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); @@ -142,8 +109,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); - printf("rShock (code units) = %g \n", rShock); - auto geom = Geometry::GetCoordinateSystem(rc); pmb->par_for( @@ -154,6 +119,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x3 = coords.Xc<3>(k, j, i); Real r = std::abs(x1); + const Real gamma = 4. / 3.; + const Real alpha0 = std::sqrt(1. - (rs / rShock)); + const Real W0 = 1. / alpha0; + const Real vr0 = -1. * std::sqrt(std::pow(W0, 2) - 1.) / W0; + const Real epsND = 0.003; // eps_bar * (W0-1), eps_bar=0.3 + const Real eta = 0.95; + Real eos_lambda[2]; if (iye > 0) { v(iye, k, j, i) = 0.5; @@ -163,38 +135,29 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // postshock - 1 if (r < rShock) { - const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); - const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); - const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - const Real T0 = temperature_from_rho_mach(eos, rho_0_Shock, target_mach, Tmin, - Tmax, vr0, eos_lambda[0]); - const Real p0 = eos.PressureFromDensityTemperature(rho_0_Shock, T0, eos_lambda); - const Real eps0 = - eos.InternalEnergyFromDensityTemperature(rho_0_Shock, T0, eos_lambda); - const Real h0 = 1. + eps0 + p0 / rho_0_Shock; const Real alphasq = 1. - (rs / r); - const Real W1 = 1.; - - Real rho1 = density_from_jump_conditions(alphasq, p0, h0, vr0, W0, rho_0_Shock, - W1, rhomin, rhomax); - const Real p1 = - (-h0 * std::pow(vr0, 2) * std::pow(W0, 3) * std::pow(rho_0_Shock, 2) + - alphasq * p0 * W1 * rho1 + - h0 * std::pow(vr0, 2) * std::pow(W0, 2) * W1 * rho_0_Shock * rho1) / - (alphasq * W1 * rho1); - const Real vr1 = (vr0 * W0 * rho_0_Shock) / (W1 * rho1); - const Real h1 = - (h0 * vr0 * std::pow(W0, 2) * rho_0_Shock) / (vr1 * std::pow(W1, 2) * rho1); - const Real eps1 = h1 - 1. - p1 / rho1; + const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1. - epsND) / W0); + const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; + const Real rho1 = (rho_0_Shock * W0 * vr0) / (vr1); + Real eps1 = (W0 - 1. + (gamma - 1) * epsND) / gamma; + + if (eps1 <= epsND) { + eps1 = 0; + } else if (eps1 > epsND && eps1 <= epsND * (eta + 1.) / eta) { + eps1 = eps1 - (eta * (eps1 - epsND)); + } else { + eps1 = eps1 - epsND; + } v(irho, k, j, i) = rho1; v(itmp, k, j, i) = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); - v(ieng, k, j, i) = rho1 * eps1; - v(iprs, k, j, i) = p1; + v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature( + rho1, v(itmp, k, j, i), eos_lambda); + v(iprs, k, j, i) = + eos.PressureFromDensityTemperature(rho1, v(itmp, k, j, i), eos_lambda); v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); @@ -216,18 +179,14 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // preshock - 0 } else { - const Real lapse0 = geom.Lapse(CellLocation::Cent, k, j, rShock); - const Real vr0 = -1. * std::sqrt(lapse0 * lapse0 - 1.); - const Real W0 = 1. / sqrt(1. - std::pow(vr0, 2)); const Real rho_0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - Real T = temperature_from_rho_mach(eos, rho_0, target_mach, Tmin, Tmax, vr0, - eos_lambda[0]); - + Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, + eos_lambda[0]); v(irho, k, j, i) = rho_0; - v(itmp, k, j, i) = T; - v(ieng, k, j, i) = - rho_0 * eos.InternalEnergyFromDensityTemperature(rho_0, T, eos_lambda); + v(itmp, k, j, i) = + eos.TemperatureFromDensityInternalEnergy(rho_0, eps0, eos_lambda); + v(ieng, k, j, i) = rho_0 * eps0; v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( @@ -248,24 +207,19 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } } - - int num_nans = std::isnan(v(ivlo, k, j, i)) + std::isnan(v(irho, k, j, i)) + - std::isnan(v(ieng, k, j, i)) + std::isnan(v(itmp, k, j, i)) + - std::isnan(v(iprs, k, j, i)) + std::isnan(v(igm1, k, j, i)); - PARTHENON_REQUIRE(num_nans == 0, "NaNs founds"); }); - fluid::PrimitiveToConserved(rc); } KOKKOS_FUNCTION -Real temperature_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real Tmin, const Real Tmax, const Real vr0, - const Real Ye) { +Real eps_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, + const Real epsmin, const Real epsmax, const Real vr0, + const Real Ye) { MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; - Real Troot = root.regula_falsi(res, Tmin, Tmax, 1.e-6 * target_mach, Tmin - 1e-10); - return Troot; + Real epsroot = + root.regula_falsi(res, epsmin, epsmax, 1.e-6 * target_mach, epsmin * 1.2); + return epsroot; } KOKKOS_FUNCTION @@ -282,14 +236,4 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]) { return (-BB - std::sqrt(discr)) / (2. * AA); } -KOKKOS_FUNCTION -Real density_from_jump_conditions(const Real alphasq, const Real p0, const Real h0, - const Real v0, const Real W0, const Real rho0, - const Real W1, const Real rhomin, const Real rhomax) { - EnergyResidual res(alphasq, p0, h0, v0, W0, rho0, W1); - root_find::RootFind root; - Real rho1root = root.regula_falsi(res, rhomin, rhomax, 1.e-10, rhomin * 5); - return rho1root; -} - } // namespace standing_accretion_shock From f7997a1aa84a824c133538bdad5b0de6ba4eb2ae Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 27 Jun 2023 16:08:31 -0600 Subject: [PATCH 39/51] fix postshock density --- src/pgen/standing_accretion_shock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 0234d35fc..b2d4506f2 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -140,7 +140,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real alphasq = 1. - (rs / r); const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1. - epsND) / W0); const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - const Real rho1 = (rho_0_Shock * W0 * vr0) / (vr1); + const Real rho1 = (rho_0 * W0 * vr0) / (vr1); Real eps1 = (W0 - 1. + (gamma - 1) * epsND) / gamma; if (eps1 <= epsND) { From 128934665f6f2ffcda4fd8896b84faccbc8ae6ac Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 27 Jun 2023 16:35:42 -0600 Subject: [PATCH 40/51] actually fix it --- src/pgen/standing_accretion_shock.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index b2d4506f2..c2d986c29 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -135,6 +135,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // postshock - 1 if (r < rShock) { + const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); const Real alphasq = 1. - (rs / r); From 51613438de962bebe72039e264502f609063bac1 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 27 Jun 2023 17:15:43 -0600 Subject: [PATCH 41/51] switch to FMKS, BL doesnt work? --- inputs/standing_accretion_shock.pin | 30 ++++++++++++++++----------- src/pgen/standing_accretion_shock.cpp | 18 ++++++++++++---- 2 files changed, 32 insertions(+), 16 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index b91d9e160..c23ac5502 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -35,23 +35,23 @@ variables = p.density, & g.n.coord file_type = hdf5 # Tabular data dump -dt = 5.e-3 # time increment between outputs +dt = 100 # time increment between outputs -nlim = 100 # cycle limit +nlim = -1 # cycle limit tlim = 5000 # time limit -integrator = rk4 # time integration algorithm -ncycle_out = 1 # interval for stdout summary info -dt_init_fact = 1.e-9 +integrator = rk2 # time integration algorithm +ncycle_out = 100 # interval for stdout summary info +dt_init_fact = 1.e-6 nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 256 # Number of zones in X1-direction -x1min = 31.2942 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 -x1max = 246.087 # maximum value of X1 | rmax | x1max => 472.5 (km) +nx1 = 1024 # Number of zones in X1-direction +x1min = 3.44199 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 +x1max = 5.43896 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag @@ -77,7 +77,7 @@ num_threads = 1 # maximum number of OMP threads field = c.c.bulk.rho method = derivative_order_1 -max_level = 3 +max_level = 1 type = IdealGas @@ -102,17 +102,23 @@ Ye = true mhd = false -enable_floors = false +enable_floors = true enable_ceilings = false -enable_flux_fixup = false -enable_c2p_fixup = false +enable_flux_fixup = true +enable_c2p_fixup = true floor_type = ConstantRhoSie rho0_floor = 1.e-3 sie0_floor = 1.e-8 +rho0_ceiling = 5.e1 +sie0_ceiling = 1.e1 + a = 0 + +derefine_poles = false + scale_free = false geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index c2d986c29..e3585e0eb 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -57,9 +57,8 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]); void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - PARTHENON_REQUIRE( - typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::BoyerLindquist), - "Problem \"standing_accretion_shock\" requires \"BoyerLindquist\" geometry!"); + PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), + "Problem \"standing_accretion_shock\" requires \"FMKS\" geometry!"); auto rc = pmb->meshblock_data.Get().get(); @@ -111,6 +110,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto geom = Geometry::GetCoordinateSystem(rc); + // set up transformation stuff + auto gpkg = pmb->packages.Get("geometry"); + bool derefine_poles = gpkg->Param("derefine_poles"); + Real h = gpkg->Param("h"); + Real xt = gpkg->Param("xt"); + Real alpha = gpkg->Param("alpha"); + Real x0 = gpkg->Param("x0"); + Real smooth = gpkg->Param("smooth"); + auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + pmb->par_for( "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { @@ -118,7 +127,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); - Real r = std::abs(x1); + // Real r = std::abs(x1); + Real r = tr.bl_radius(x1); const Real gamma = 4. / 3.; const Real alpha0 = std::sqrt(1. - (rs / rShock)); const Real W0 = 1. / alpha0; From 2adc4252c6aa57a6e43df9918d26a57510d042bb Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 28 Jun 2023 11:59:46 -0600 Subject: [PATCH 42/51] set min eps and T --- src/pgen/standing_accretion_shock.cpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index e3585e0eb..19b0475f8 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -98,15 +98,18 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - - const Real epsmin = 1.e-9; - const Real epsmax = 1.e6; + auto epsmin = pmb->packages.Get("eos")->Param("sie_min"); + auto epsmax = pmb->packages.Get("eos")->Param("sie_max"); + auto Cv = pmb->packages.Get("eos")->Param("Cv"); + auto Tmin = pmb->packages.Get("eos")->Param("T_min"); + auto Tmax = pmb->packages.Get("eos")->Param("T_max"); Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); Real MPNS = 1.3 * solar_mass; Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + i epsmin *= (-1.); auto geom = Geometry::GetCoordinateSystem(rc); @@ -155,7 +158,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real eps1 = (W0 - 1. + (gamma - 1) * epsND) / gamma; if (eps1 <= epsND) { - eps1 = 0; + eps1 = epsmin; } else if (eps1 > epsND && eps1 <= epsND * (eta + 1.) / eta) { eps1 = eps1 - (eta * (eps1 - epsND)); } else { @@ -195,8 +198,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, eos_lambda[0]); v(irho, k, j, i) = rho_0; - v(itmp, k, j, i) = - eos.TemperatureFromDensityInternalEnergy(rho_0, eps0, eos_lambda); + v(itmp, k, j, i) = std::max(Tmin, eps0 / Cv); // eos.TemperatureFromDensityInternalEnergy(rho_0, + // eps0, eos_lambda); v(ieng, k, j, i) = rho_0 * eps0; v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); From 8ec33b5487dc1ddce622a87228ec0a2a47c8dfdc Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 28 Jun 2023 12:06:44 -0600 Subject: [PATCH 43/51] i,i,i --- src/pgen/standing_accretion_shock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 19b0475f8..c45e63816 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -109,7 +109,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real MPNS = 1.3 * solar_mass; Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); - i epsmin *= (-1.); + epsmin *= (-1.); auto geom = Geometry::GetCoordinateSystem(rc); From d2172aeb755a2ed3b5829b80ee91809b091ae107 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 29 Jun 2023 09:26:52 -0600 Subject: [PATCH 44/51] formattttttttttttttttttt --- src/pgen/standing_accretion_shock.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index c45e63816..e3d202bb1 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -98,8 +98,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - auto epsmin = pmb->packages.Get("eos")->Param("sie_min"); - auto epsmax = pmb->packages.Get("eos")->Param("sie_max"); + const Real epsmin = 1.e-10; + const Real epsmax = 1.e10; auto Cv = pmb->packages.Get("eos")->Param("Cv"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); @@ -109,7 +109,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real MPNS = 1.3 * solar_mass; Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); - epsmin *= (-1.); + epsmin = std::abs(epsmin); auto geom = Geometry::GetCoordinateSystem(rc); @@ -198,8 +198,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, eos_lambda[0]); v(irho, k, j, i) = rho_0; - v(itmp, k, j, i) = std::max(Tmin, eps0 / Cv); // eos.TemperatureFromDensityInternalEnergy(rho_0, - // eps0, eos_lambda); + v(itmp, k, j, i) = std::max( + Tmin, eps0 / Cv); // eos.TemperatureFromDensityInternalEnergy(rho_0, + // eps0, eos_lambda); v(ieng, k, j, i) = rho_0 * eps0; v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); From 3906a331a5078b259d139a60723ca190d4335ad2 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 5 Jul 2023 12:48:23 -0600 Subject: [PATCH 45/51] debugging --- inputs/standing_accretion_shock.pin | 34 +++++------- src/pgen/standing_accretion_shock.cpp | 77 +++++++++++---------------- 2 files changed, 42 insertions(+), 69 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index c23ac5502..ffaed0fa1 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -49,9 +49,9 @@ nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 1024 # Number of zones in X1-direction -x1min = 3.44199 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 -x1max = 5.43896 # maximum value of X1 | rmax | x1max => 472.5 (km) +nx1 = 512 # Number of zones in X1-direction +x1min = 31.2492 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 +x1max = 230.202 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag @@ -79,10 +79,15 @@ field = c.c.bulk.rho method = derivative_order_1 max_level = 1 +# +#type = IdealGas +#Gamma = 1.33 +#Cv = 3.0 + -type = IdealGas -Gamma = 1.3333333333333333 -Cv = 3.0 +type = StellarCollapse +filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 +use_sp5 = false hydro = true @@ -92,7 +97,7 @@ rad = false xorder = 2 -cfl = 0.5 +cfl = 0.6 riemann = llf recon = weno5 c2p_method = robust @@ -101,24 +106,9 @@ c2p_tol = 1.e-8 Ye = true mhd = false - -enable_floors = true -enable_ceilings = false -enable_flux_fixup = true -enable_c2p_fixup = true -floor_type = ConstantRhoSie -rho0_floor = 1.e-3 -sie0_floor = 1.e-8 -rho0_ceiling = 5.e1 -sie0_ceiling = 1.e1 - - a = 0 - -derefine_poles = false - scale_free = false geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index e3d202bb1..899a8d588 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -14,7 +14,6 @@ // Parthenon #include -#include "geometry/boyer_lindquist.hpp" #include "pgen/pgen.hpp" #include "phoebus_utils/root_find.hpp" #include "phoebus_utils/unit_conversions.hpp" @@ -57,10 +56,10 @@ Real ucon_norm(Real ucon[4], Real gcov[4][4]); void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::FMKS), - "Problem \"standing_accretion_shock\" requires \"FMKS\" geometry!"); + PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::SphericalKerrSchild), + "Problem \"standing_accretion_shock\" requires \"SphericalKerrSchild\" geometry!"); - auto rc = pmb->meshblock_data.Get().get(); + auto &rc = pmb->meshblock_data.Get(); PackIndexMap imap; auto v = rc->PackVariables( @@ -81,8 +80,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const std::string eos_type = pin->GetString("eos", "type"); PARTHENON_REQUIRE_THROWS( - eos_type == "IdealGas", - "Standing Accretion Shock setup only works with Ideal Gas EOS"); + eos_type == "IdealGas" || eos_type == "StellarCollapse", + "Standing Accretion Shock setup only works with Ideal Gas or Stellar Collapse EOS"); Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); @@ -98,31 +97,21 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto eos = pmb->packages.Get("eos")->Param("d.EOS"); const Real a = pin->GetReal("geometry", "a"); auto bl = Geometry::BoyerLindquist(a); - const Real epsmin = 1.e-10; - const Real epsmax = 1.e10; - auto Cv = pmb->packages.Get("eos")->Param("Cv"); + auto epsmin = pmb->packages.Get("eos")->Param("sie_min"); + auto epsmax = pmb->packages.Get("eos")->Param("sie_max"); + //auto Cv = pmb->packages.Get("eos")->Param("Cv"); auto Tmin = pmb->packages.Get("eos")->Param("T_min"); auto Tmax = pmb->packages.Get("eos")->Param("T_max"); + printf("Tmin, Tmax, epsmin, epsmax = %g %g %g %g\n", Tmin,Tmax, epsmin, epsmax); + Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); Real MPNS = 1.3 * solar_mass; - Real rs = - (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); - epsmin = std::abs(epsmin); - - auto geom = Geometry::GetCoordinateSystem(rc); - - // set up transformation stuff - auto gpkg = pmb->packages.Get("geometry"); - bool derefine_poles = gpkg->Param("derefine_poles"); - Real h = gpkg->Param("h"); - Real xt = gpkg->Param("xt"); - Real alpha = gpkg->Param("alpha"); - Real x0 = gpkg->Param("x0"); - Real smooth = gpkg->Param("smooth"); - auto tr = Geometry::McKinneyGammieRyan(derefine_poles, h, xt, alpha, x0, smooth); + Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + auto geom = Geometry::GetCoordinateSystem(rc.get()); + printf("Rs, rmin (code units) = %g %g \n", rs, std::abs(coords.Xc<1>(1, 1, 1))); pmb->par_for( "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { @@ -130,8 +119,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); - // Real r = std::abs(x1); - Real r = tr.bl_radius(x1); + Real r = std::abs(x1); const Real gamma = 4. / 3.; const Real alpha0 = std::sqrt(1. - (rs / rShock)); const Real W0 = 1. / alpha0; @@ -139,7 +127,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real epsND = 0.003; // eps_bar * (W0-1), eps_bar=0.3 const Real eta = 0.95; - Real eos_lambda[2]; + Real eos_lambda[2] = {0.}; if (iye > 0) { v(iye, k, j, i) = 0.5; eos_lambda[0] = v(iye, k, j, i); @@ -149,8 +137,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (r < rShock) { const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - const Real rho_0_Shock = - Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); const Real alphasq = 1. - (rs / r); const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1. - epsND) / W0); const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; @@ -165,16 +152,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eps1 = eps1 - epsND; } + Real T1 = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); + v(irho, k, j, i) = rho1; - v(itmp, k, j, i) = - eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); - v(ieng, k, j, i) = rho1 * eos.InternalEnergyFromDensityTemperature( - rho1, v(itmp, k, j, i), eos_lambda); - v(iprs, k, j, i) = - eos.PressureFromDensityTemperature(rho1, v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); + v(itmp, k, j, i) = T1; + v(ieng, k, j, i) = rho1 * eps1; + v(iprs, k, j, i) = eos.PressureFromDensityTemperature(rho1, v(itmp, k, j, i), eos_lambda); + v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); Real ucon[] = {0.0, vr1, 0.0, 0.0}; Real gcov[4][4]; @@ -193,14 +177,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // preshock - 0 } else { - const Real rho_0 = - Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, - eos_lambda[0]); + const Real rho_0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); + const Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, eos_lambda[0]); + Real T0 = eos.TemperatureFromDensityInternalEnergy(rho_0,eps0, eos_lambda); + v(irho, k, j, i) = rho_0; - v(itmp, k, j, i) = std::max( - Tmin, eps0 / Cv); // eos.TemperatureFromDensityInternalEnergy(rho_0, - // eps0, eos_lambda); + v(itmp, k, j, i) = T0; v(ieng, k, j, i) = rho_0 * eps0; v(iprs, k, j, i) = eos.PressureFromDensityTemperature( v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); @@ -222,8 +204,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } } + printf("r, rho, T, eps, W0, eos_lambda[0], eos_lambda[1] = %g %g %g %g %g %g %g \n", r, v(irho, k, j, i), v(itmp, k, j, i), v(ieng, k, j, i) / v(irho, k, j, i), W0, eos_lambda[0], eos_lambda[1]); }); - fluid::PrimitiveToConserved(rc); + fluid::PrimitiveToConserved(rc.get()); } KOKKOS_FUNCTION @@ -233,7 +216,7 @@ Real eps_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, MachResidual res(eos, rho, vr0, target_mach, Ye); root_find::RootFind root; Real epsroot = - root.regula_falsi(res, epsmin, epsmax, 1.e-6 * target_mach, epsmin * 1.2); + root.regula_falsi(res, epsmin, epsmax, 1.e-6 * target_mach, std::max(epsmin,1e-10)); return epsroot; } From 6412ce1e82a76084f54f9f4efc910ac0f59269e2 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 18 Jul 2023 16:02:13 -0600 Subject: [PATCH 46/51] reads in ics for sasw --- inputs/standing_accretion_shock.pin | 44 +++-- src/pgen/standing_accretion_shock.cpp | 237 ++++++++++++-------------- 2 files changed, 134 insertions(+), 147 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index ffaed0fa1..44e99d275 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -35,11 +35,11 @@ variables = p.density, & g.n.coord file_type = hdf5 # Tabular data dump -dt = 100 # time increment between outputs +dt = 10 # time increment between outputs nlim = -1 # cycle limit -tlim = 5000 # time limit +tlim = 2000 # time limit integrator = rk2 # time integration algorithm ncycle_out = 100 # interval for stdout summary info dt_init_fact = 1.e-6 @@ -49,9 +49,9 @@ nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 512 # Number of zones in X1-direction -x1min = 31.2492 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 -x1max = 230.202 # maximum value of X1 | rmax | x1max => 472.5 (km) +nx1 = 1024 # Number of zones in X1-direction +x1min = 30 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 +x1max = 230 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag ox1_bc = user # Outer-X1 boundary condition flag @@ -79,15 +79,28 @@ field = c.c.bulk.rho method = derivative_order_1 max_level = 1 -# -#type = IdealGas -#Gamma = 1.33 -#Cv = 3.0 - -type = StellarCollapse -filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 -use_sp5 = false +type = IdealGas +Gamma = 1.33 +Cv = 3.0 + +# +#type = StellarCollapse +#filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 +#use_sp5 = false + + +enable_floors = true +enable_mhd_floors = false +enable_ceilings = true +enable_flux_fixup = true +enable_c2p_fixup = true +floor_type = ConstantRhoSie +rho0_floor = 1.e-12 +sie0_floor = 1.e-6 +ceiling_type = ConstantGamSie +gam0_ceiling = 100.0 +sie_ceiling = 1e-2 hydro = true @@ -111,10 +124,11 @@ a = 0 scale_free = false -geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km -fluid_density_cgs = 1.e9 # gives mass = 7.07e24 g +geom_mass_msun = 1.3 # canonical M=1.3Msun NS sets our length scale of L = G M / c^2 = 1.93 km +fluid_density_cgs = 3.599e17 # 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km target_mach = -100 # target Mach number for upstream flow in preshock region +model_filename = 1d_initial_conditions_radice.txt diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index 899a8d588..a1a41c719 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -1,4 +1,5 @@ -// © 2021. Triad National Security, LLC. All rights reserved. This +// parthenon::ParArray2D © 2021. Triad National Security, LLC. All rights reserved. +// This // program was produced under U.S. Government contract // 89233218CNA000001 for Los Alamos National Laboratory (LANL), which // is operated by Triad National Security, LLC for the U.S. @@ -14,50 +15,43 @@ // Parthenon #include +#include "monopole_gr/monopole_gr_base.hpp" #include "pgen/pgen.hpp" #include "phoebus_utils/root_find.hpp" #include "phoebus_utils/unit_conversions.hpp" #include "utils/error_checking.hpp" #include +// Spiner +#include + namespace standing_accretion_shock { parthenon::constants::PhysicalConstants pc; using Microphysics::EOS::EOS; -class MachResidual { - public: - KOKKOS_FUNCTION - MachResidual(const EOS &eos, const Real rho, const Real vr0, const Real target_mach, - const Real Ye) - : eos_(eos), rho_(rho), vr0_(vr0), target_mach_(target_mach) { - lambda_[0] = Ye; - } - KOKKOS_INLINE_FUNCTION - Real operator()(const Real eps) { - const Real gamma = 4. / 3.; - Real cs = std::sqrt(gamma * (gamma - 1.) * eps); - Real mach = vr0_ / cs; - return mach - target_mach_; - } +constexpr int NSASW = 4; +constexpr int R = 0; +constexpr int RHO = 1; +constexpr int VR = 2; +constexpr int EPS = 3; - private: - const EOS &eos_; - Real rho_, vr0_, target_mach_; - Real lambda_[2]; -}; +using State_t = parthenon::ParArray2D; +using State_host_t = typename parthenon::ParArray2D::HostMirror; +using Radius = Spiner::RegularGrid1D; + +template +void Get1DProfileData(const std::string model_filename, const int num_zones, + const int num_comments, const int num_vars, ParArray2D &state_raw); -KOKKOS_FUNCTION -Real eps_from_rho_mach(const Microphysics::EOS::EOS &eos, const Real rho, - const Real target_mach, const Real epsmin, const Real epsmax, - const Real vr0, const Real Ye); KOKKOS_FUNCTION Real ucon_norm(Real ucon[4], Real gcov[4][4]); void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { - PARTHENON_REQUIRE(typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::SphericalKerrSchild), - "Problem \"standing_accretion_shock\" requires \"SphericalKerrSchild\" geometry!"); + PARTHENON_REQUIRE( + typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::SphericalKerrSchild), + "Problem \"standing_accretion_shock\" requires \"SphericalKerrSchild\" geometry!"); auto &rc = pmb->meshblock_data.Get(); @@ -83,10 +77,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos_type == "IdealGas" || eos_type == "StellarCollapse", "Standing Accretion Shock setup only works with Ideal Gas or Stellar Collapse EOS"); - Real Mdot = pin->GetOrAddReal("standing_accretion_shock", "Mdot", 0.2); - Real rShock = pin->GetOrAddReal("standing_accretion_shock", "rShock", 200); - Real target_mach = pin->GetOrAddReal("standing_accretion_shock", "target_mach", -100); - IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); @@ -96,22 +86,36 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { pmb->packages.Get("phoebus")->Param("unit_conv"); auto eos = pmb->packages.Get("eos")->Param("d.EOS"); const Real a = pin->GetReal("geometry", "a"); - auto bl = Geometry::BoyerLindquist(a); - auto epsmin = pmb->packages.Get("eos")->Param("sie_min"); - auto epsmax = pmb->packages.Get("eos")->Param("sie_max"); - //auto Cv = pmb->packages.Get("eos")->Param("Cv"); - auto Tmin = pmb->packages.Get("eos")->Param("T_min"); - auto Tmax = pmb->packages.Get("eos")->Param("T_max"); + auto geom = Geometry::GetCoordinateSystem(rc.get()); - printf("Tmin, Tmax, epsmin, epsmax = %g %g %g %g\n", Tmin,Tmax, epsmin, epsmax); + // info about model file + std::string model_filename = + pin->GetOrAddString("standing_accretion_shock", "model_filename", "file.txt"); + const int num_zones = 1000 - 1; + const int num_comments = 6; + const int num_vars = 4; - Mdot *= ((solar_mass * unit_conv.GetMassCGSToCode()) / unit_conv.GetTimeCGSToCode()); - rShock *= (1.e5 * unit_conv.GetLengthCGSToCode()); - Real MPNS = 1.3 * solar_mass; - Real rs = (2. * pc.g_newt * MPNS / (std::pow(pc.c, 2))) * unit_conv.GetLengthCGSToCode(); + // Kokkos views for loading in data + State_t state_raw_d("state raw", NSASW, num_zones); + State_host_t state_raw_h = Kokkos::create_mirror_view(state_raw_d); + + Get1DProfileData(model_filename, num_zones, num_comments, num_vars, state_raw_h); + + Kokkos::deep_copy(state_raw_d, state_raw_h); + + // auto base = params.Get("state_raw_d"); + + const Real rhomin = 1.e-6; + const Real epsmin = 1e-16; + + const Real rin = state_raw_h(R, 0); + const Real rout = state_raw_h(R, num_zones - 1); + + printf("rin = %.14e (code units)\n", rin); + printf("rout = %.14e (code units)\n", rout); + + Radius raw_radius(rin, rout, num_zones); - auto geom = Geometry::GetCoordinateSystem(rc.get()); - printf("Rs, rmin (code units) = %g %g \n", rs, std::abs(coords.Xc<1>(1, 1, 1))); pmb->par_for( "Phoebus::ProblemGenerator::StandingAccrectionShock", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { @@ -120,104 +124,73 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { const Real x3 = coords.Xc<3>(k, j, i); Real r = std::abs(x1); - const Real gamma = 4. / 3.; - const Real alpha0 = std::sqrt(1. - (rs / rShock)); - const Real W0 = 1. / alpha0; - const Real vr0 = -1. * std::sqrt(std::pow(W0, 2) - 1.) / W0; - const Real epsND = 0.003; // eps_bar * (W0-1), eps_bar=0.3 - const Real eta = 0.95; - - Real eos_lambda[2] = {0.}; + Real eos_lambda[2] = {0.}; + if (iye > 0) { v(iye, k, j, i) = 0.5; eos_lambda[0] = v(iye, k, j, i); } - // postshock - 1 - if (r < rShock) { - - const Real rho_0 = Mdot / (4. * M_PI * std::pow(r, 2) * W0 * std::abs(vr0)); - const Real rho_0_Shock = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - const Real alphasq = 1. - (rs / r); - const Real psi = alphasq * ((gamma - 1.) / gamma) * ((W0 - 1. - epsND) / W0); - const Real vr1 = (vr0 + std::sqrt(vr0 * vr0 - 4. * psi)) / 2.; - const Real rho1 = (rho_0 * W0 * vr0) / (vr1); - Real eps1 = (W0 - 1. + (gamma - 1) * epsND) / gamma; - - if (eps1 <= epsND) { - eps1 = epsmin; - } else if (eps1 > epsND && eps1 <= epsND * (eta + 1.) / eta) { - eps1 = eps1 - (eta * (eps1 - epsND)); - } else { - eps1 = eps1 - epsND; - } - - Real T1 = eos.TemperatureFromDensityInternalEnergy(rho1, eps1, eos_lambda); - - v(irho, k, j, i) = rho1; - v(itmp, k, j, i) = T1; - v(ieng, k, j, i) = rho1 * eps1; - v(iprs, k, j, i) = eos.PressureFromDensityTemperature(rho1, v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature(v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / v(iprs, k, j, i); - - Real ucon[] = {0.0, vr1, 0.0, 0.0}; - Real gcov[4][4]; - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); - ucon[0] = ucon_norm(ucon, gcov); - - const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); - Real beta[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real W = lapse * ucon[0]; - - for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; - } - - // preshock - 0 - } else { - - const Real rho_0 = Mdot / (4. * M_PI * std::pow(rShock, 2) * W0 * std::abs(vr0)); - const Real eps0 = eps_from_rho_mach(eos, rho_0, target_mach, epsmin, epsmax, vr0, eos_lambda[0]); - Real T0 = eos.TemperatureFromDensityInternalEnergy(rho_0,eps0, eos_lambda); - - v(irho, k, j, i) = rho_0; - v(itmp, k, j, i) = T0; - v(ieng, k, j, i) = rho_0 * eps0; - v(iprs, k, j, i) = eos.PressureFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda); - v(igm1, k, j, i) = eos.BulkModulusFromDensityTemperature( - v(irho, k, j, i), v(itmp, k, j, i), eos_lambda) / - v(iprs, k, j, i); - - Real ucon[] = {0.0, vr0, 0.0, 0.0}; - Real gcov[4][4]; - geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); - ucon[0] = ucon_norm(ucon, gcov); - - const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); - Real beta[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real W = lapse * ucon[0]; - - for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; - } + Real rho = + std::max(rhomin, MonopoleGR::Interpolate(r, state_raw_d, raw_radius, RHO)); + Real eps = + std::max(epsmin, MonopoleGR::Interpolate(r, state_raw_d, raw_radius, EPS)); + Real vr = std::max(-1., MonopoleGR::Interpolate(r, state_raw_d, raw_radius, VR)); + Real T = eos.TemperatureFromDensityInternalEnergy(rho, eps, eos_lambda); + + v(irho, k, j, i) = rho; + v(itmp, k, j, i) = T; + v(ieng, k, j, i) = rho * eps; + v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy(rho, eps, eos_lambda); + v(igm1, k, j, i) = + eos.BulkModulusFromDensityTemperature(rho, T, eos_lambda) / v(iprs, k, j, i); + + Real ucon[] = {0.0, vr, 0.0, 0.0}; + Real gcov[4][4]; + geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); + ucon[0] = ucon_norm(ucon, gcov); + + const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); + Real beta[3]; + geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); + Real W = lapse * ucon[0]; + + for (int d = 0; d < 3; d++) { + v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; } - printf("r, rho, T, eps, W0, eos_lambda[0], eos_lambda[1] = %g %g %g %g %g %g %g \n", r, v(irho, k, j, i), v(itmp, k, j, i), v(ieng, k, j, i) / v(irho, k, j, i), W0, eos_lambda[0], eos_lambda[1]); }); fluid::PrimitiveToConserved(rc.get()); } -KOKKOS_FUNCTION -Real eps_from_rho_mach(const EOS &eos, const Real rho, const Real target_mach, - const Real epsmin, const Real epsmax, const Real vr0, - const Real Ye) { - MachResidual res(eos, rho, vr0, target_mach, Ye); - root_find::RootFind root; - Real epsroot = - root.regula_falsi(res, epsmin, epsmax, 1.e-6 * target_mach, std::max(epsmin,1e-10)); - return epsroot; +template +void Get1DProfileData(const std::string model_filename, const int num_zones, + const int num_comments, const int num_vars, ParArray2D &state_raw) { + + std::ifstream inputfile(model_filename); + + Real val = 0; + std::string line; + int line_num = 0; + + // skip comment lines + while (line_num < num_comments) { + getline(inputfile, line); + line_num++; + } + + // read file into model_1d + for (int i = 0; i < num_zones; i++) // number of zones + { + for (int j = 0; j < num_vars; j++) // number of vars + { + inputfile >> val; + state_raw(j, i) = val; + } + } + + printf("Read 1D profile into state array.\n"); + + inputfile.close(); } KOKKOS_FUNCTION From 0553e9a74f4125c99999e444676b6d9d8e79aebd Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Wed, 19 Jul 2023 12:13:58 -0600 Subject: [PATCH 47/51] fix line count, format, cleanup --- src/pgen/standing_accretion_shock.cpp | 81 ++++++++++++++++++++++----- 1 file changed, 66 insertions(+), 15 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index a1a41c719..a03c49785 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -30,16 +30,22 @@ namespace standing_accretion_shock { parthenon::constants::PhysicalConstants pc; using Microphysics::EOS::EOS; -constexpr int NSASW = 4; +constexpr int NSASW = 8; constexpr int R = 0; constexpr int RHO = 1; constexpr int VR = 2; constexpr int EPS = 3; +constexpr int PRES = 4; +constexpr int HEAT = 5; +constexpr int S = 6; +constexpr int OMEGA = 7; using State_t = parthenon::ParArray2D; using State_host_t = typename parthenon::ParArray2D::HostMirror; using Radius = Spiner::RegularGrid1D; +std::pair Get1DProfileNumZones(const std::string model_filename); + template void Get1DProfileData(const std::string model_filename, const int num_zones, const int num_comments, const int num_vars, ParArray2D &state_raw); @@ -91,25 +97,27 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { // info about model file std::string model_filename = pin->GetOrAddString("standing_accretion_shock", "model_filename", "file.txt"); - const int num_zones = 1000 - 1; - const int num_comments = 6; - const int num_vars = 4; + + std::pair p = Get1DProfileNumZones(model_filename); + const int num_zones = p.first - 1; + const int num_comments = p.second; + + printf("1D model read in with %d number of zones. \n", num_zones); + printf("1D model read in with %d number of comments. \n", num_comments); // Kokkos views for loading in data State_t state_raw_d("state raw", NSASW, num_zones); State_host_t state_raw_h = Kokkos::create_mirror_view(state_raw_d); - Get1DProfileData(model_filename, num_zones, num_comments, num_vars, state_raw_h); + Get1DProfileData(model_filename, num_zones, num_comments, NSASW, state_raw_h); Kokkos::deep_copy(state_raw_d, state_raw_h); - // auto base = params.Get("state_raw_d"); - - const Real rhomin = 1.e-6; - const Real epsmin = 1e-16; + const Real rhomin = 1.e-16; + const Real epsmin = 1.e-16; - const Real rin = state_raw_h(R, 0); - const Real rout = state_raw_h(R, num_zones - 1); + Real rin = state_raw_h(R, 0); + Real rout = state_raw_h(R, num_zones - 1); printf("rin = %.14e (code units)\n", rin); printf("rout = %.14e (code units)\n", rout); @@ -136,14 +144,13 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { Real eps = std::max(epsmin, MonopoleGR::Interpolate(r, state_raw_d, raw_radius, EPS)); Real vr = std::max(-1., MonopoleGR::Interpolate(r, state_raw_d, raw_radius, VR)); - Real T = eos.TemperatureFromDensityInternalEnergy(rho, eps, eos_lambda); v(irho, k, j, i) = rho; - v(itmp, k, j, i) = T; v(ieng, k, j, i) = rho * eps; v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy(rho, eps, eos_lambda); v(igm1, k, j, i) = - eos.BulkModulusFromDensityTemperature(rho, T, eos_lambda) / v(iprs, k, j, i); + eos.BulkModulusFromDensityInternalEnergy(rho, eps, eos_lambda) / + v(iprs, k, j, i); Real ucon[] = {0.0, vr, 0.0, 0.0}; Real gcov[4][4]; @@ -162,6 +169,50 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { fluid::PrimitiveToConserved(rc.get()); } +std::pair Get1DProfileNumZones(const std::string model_filename) { + + // open file + std::ifstream inputfile(model_filename); + const std::string whitespace(" \t\n\r"); + + // error check + if (!inputfile.is_open()) { + std::cout << model_filename << " not found :( \n."; + exit(-1); + } + + int nz = 0; + int nc = 0; + std::string line; + std::string comment1("#"); + std::string comment2("//"); + + // get number of zones from 1d file + while (!inputfile.eof()) { + + getline(inputfile, line); + + std::size_t first_nonws = line.find_first_not_of(whitespace); + + // skip empty lines + if (first_nonws == std::string::npos) { + continue; + } + + // skip comments + if (line.find(comment1) == first_nonws || line.find(comment2) == first_nonws) { + nc++; + continue; + } + + nz++; + } + + inputfile.close(); + + return std::make_pair(nz + 1, nc); +} + template void Get1DProfileData(const std::string model_filename, const int num_zones, const int num_comments, const int num_vars, ParArray2D &state_raw) { @@ -188,7 +239,7 @@ void Get1DProfileData(const std::string model_filename, const int num_zones, } } - printf("Read 1D profile into state array.\n"); + printf("Read 1D profile into state array.\n"); inputfile.close(); } From 4a41d06dcea2bb15af857993d8b8a4be887f2f55 Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Thu, 20 Jul 2023 12:24:10 -0600 Subject: [PATCH 48/51] format, fix velocity --- src/pgen/standing_accretion_shock.cpp | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp index a03c49785..7dcff2ab7 100644 --- a/src/pgen/standing_accretion_shock.cpp +++ b/src/pgen/standing_accretion_shock.cpp @@ -35,10 +35,6 @@ constexpr int R = 0; constexpr int RHO = 1; constexpr int VR = 2; constexpr int EPS = 3; -constexpr int PRES = 4; -constexpr int HEAT = 5; -constexpr int S = 6; -constexpr int OMEGA = 7; using State_t = parthenon::ParArray2D; using State_host_t = typename parthenon::ParArray2D::HostMirror; @@ -152,18 +148,20 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { eos.BulkModulusFromDensityInternalEnergy(rho, eps, eos_lambda) / v(iprs, k, j, i); - Real ucon[] = {0.0, vr, 0.0, 0.0}; + // three-velocity + const Real vcon[3] = {vr, 0.0, 0.0}; + Real vsq = 0.; + SPACELOOP2(ii, jj) { vsq += vcon[ii] * vcon[jj]; } + const Real W = 1. / sqrt(1. - vsq); + + // four-velocity + Real ucon[] = {0.0, W * vcon[0], 0.0, 0.0}; Real gcov[4][4]; geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov); ucon[0] = ucon_norm(ucon, gcov); - const Real lapse = geom.Lapse(CellLocation::Cent, k, j, i); - Real beta[3]; - geom.ContravariantShift(CellLocation::Cent, k, j, i, beta); - Real W = lapse * ucon[0]; - for (int d = 0; d < 3; d++) { - v(ivlo + d, k, j, i) = ucon[d + 1] + W * beta[d] / lapse; + v(ivlo + d, k, j, i) = ucon[d + 1]; } }); fluid::PrimitiveToConserved(rc.get()); From b26325c2a63448d0e36e84e417697cf167007ac0 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Fri, 6 Oct 2023 15:57:58 -0400 Subject: [PATCH 49/51] cartesian coordinate transformation added in cooling_function for lightbulb --- src/pgen/progenitor.cpp | 9 ++++----- src/radiation/cooling_function.cpp | 22 +++++++++++++++++++++- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/pgen/progenitor.cpp b/src/pgen/progenitor.cpp index 0c6cda5d5..d5a24a127 100644 --- a/src/pgen/progenitor.cpp +++ b/src/pgen/progenitor.cpp @@ -115,17 +115,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (is_monopole_sph) { r = std::abs(x1); } else { // Cartesian - r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3); - transform(x1, x2, x3, C, c2s, s2c); - } + r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3); + transform(x1, x2, x3, C, c2s, s2c); + } Real lambda[2]; if (iye > 0) { v(iye, k, j, i) = Ye_dev.interpToReal(r); lambda[0] = v(iye, k, j, i); } - - const Real u = phoebus::energy_from_rho_P(eos, mass_density_dev.interpToReal(r), + const Real u = phoebus::energy_from_rho_P(eos, mass_density_dev.interpToReal(r), pressure_dev.interpToReal(r), emin, emax, lambda[0]); const Real sie = u / mass_density_dev.interpToReal(r); diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp index 642618e82..02e4bb7ba 100644 --- a/src/radiation/cooling_function.cpp +++ b/src/radiation/cooling_function.cpp @@ -165,7 +165,15 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshBlockData *rc, const doub }); // Light Bulb with Liebendorfer model + const bool is_monopole_cart = + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); + const bool is_monopole_sph = + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph)); + auto &coords = pmb->coords; + using Transformation_t = Geometry::SphericalToCartesian; + auto const &geom_pkg = pmb->packages.Get("geometry"); + auto transform = Geometry::GetTransformation(geom_pkg.get()); const bool do_liebendorfer = rad->Param("do_liebendorfer"); const bool do_lightbulb = rad->Param("do_lightbulb"); if (do_lightbulb) { @@ -180,7 +188,19 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshBlockData *rc, const doub DEFAULT_LOOP_PATTERN, "CoolingFunctionCalculateFourForce", DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { - const Real r = std::abs(coords.Xc<1>(k, j, i)); // TODO(MG) coord transform game + const Real x1 = coords.Xc<1>(k, j, i); + const Real x2 = coords.Xc<2>(k, j, i); + const Real x3 = coords.Xc<3>(k, j, i); + Real Cart[3]; + Real s2c[3][3], c2s[3][3]; + Real r; + + if (is_monopole_sph){ + r = std::abs(x1); + } else { + r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3); + transform(x1, x2, x3, Cart, c2s, s2c); + } const Real rho = v(prho, k, j, i) * unit_conv.GetMassDensityCodeToCGS(); // Density in CGS const Real cdensity = v(crho, k, j, i); // conserved density From e08387b1bb8a4f1804785d0379c7712bd12b794e Mon Sep 17 00:00:00 2001 From: carlnotsagan Date: Tue, 7 Nov 2023 13:24:04 -0700 Subject: [PATCH 50/51] try to turn on lightbulb for SASW using radice ICs --- inputs/standing_accretion_shock.pin | 47 ++++++++++++++++++----------- 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin index 44e99d275..b25e6c3e0 100644 --- a/inputs/standing_accretion_shock.pin +++ b/inputs/standing_accretion_shock.pin @@ -32,24 +32,27 @@ variables = p.density, & temperature,& p.ye, & g.c.coord, & - g.n.coord + g.n.coord, & + flux_divergence, & + src_terms + file_type = hdf5 # Tabular data dump -dt = 10 # time increment between outputs +dt = 100 # time increment between outputs nlim = -1 # cycle limit -tlim = 2000 # time limit +tlim = 10000 # time limit integrator = rk2 # time integration algorithm ncycle_out = 100 # interval for stdout summary info -dt_init_fact = 1.e-6 +dt_init_fact = 1.e-8 nghost = 4 #refinement = adaptive #numlevel = 3 -nx1 = 1024 # Number of zones in X1-direction +nx1 = 512 # Number of zones in X1-direction x1min = 30 # minimum value of X1 | rPNS | x1min = ln R (code units), R = 60 (km) => R = 31.2492 (code units) => x1min = 3.44199 x1max = 230 # maximum value of X1 | rmax | x1max => 472.5 (km) ix1_bc = user # Inner-X1 boundary condition flag @@ -79,21 +82,21 @@ field = c.c.bulk.rho method = derivative_order_1 max_level = 1 - -type = IdealGas -Gamma = 1.33 -Cv = 3.0 - # -#type = StellarCollapse -#filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 -#use_sp5 = false +#type = IdealGas +#Gamma = 1.33 +#Cv = 3.0 + + +type = StellarCollapse +filename = /usr/projects/w22_stars/EoS/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5 +use_sp5 = false -enable_floors = true +enable_floors = false enable_mhd_floors = false -enable_ceilings = true -enable_flux_fixup = true +enable_ceilings = false +enable_flux_fixup = false enable_c2p_fixup = true floor_type = ConstantRhoSie rho0_floor = 1.e-12 @@ -106,7 +109,7 @@ sie_ceiling = 1e-2 hydro = true he = false 3t = false -rad = false +rad = true xorder = 2 @@ -131,4 +134,12 @@ fluid_density_cgs = 3.599e17 # 1.e9 # gives mass = 7.07e24 g Mdot = 0.2 # Msun / sec rShock = 200 # km target_mach = -100 # target Mach number for upstream flow in preshock region -model_filename = 1d_initial_conditions_radice.txt +model_filename = 1d_initial_conditions_radice_10k.txt + + +method=cooling_function +do_liebendorfer=false +do_lightbulb=true +lum=2.5 +do_gain_reducer=false +always_gain=false From 2d825c7d3647ef8d0d62ed2da9ef248b8c54bbf3 Mon Sep 17 00:00:00 2001 From: mari2895 Date: Wed, 8 Nov 2023 17:50:12 -0500 Subject: [PATCH 51/51] formatting applied --- src/pgen/progenitor.cpp | 8 ++++---- src/radiation/cooling_function.cpp | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/pgen/progenitor.cpp b/src/pgen/progenitor.cpp index d5a24a127..b6409e008 100644 --- a/src/pgen/progenitor.cpp +++ b/src/pgen/progenitor.cpp @@ -115,16 +115,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (is_monopole_sph) { r = std::abs(x1); } else { // Cartesian - r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3); - transform(x1, x2, x3, C, c2s, s2c); - } + r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3); + transform(x1, x2, x3, C, c2s, s2c); + } Real lambda[2]; if (iye > 0) { v(iye, k, j, i) = Ye_dev.interpToReal(r); lambda[0] = v(iye, k, j, i); } - const Real u = phoebus::energy_from_rho_P(eos, mass_density_dev.interpToReal(r), + const Real u = phoebus::energy_from_rho_P(eos, mass_density_dev.interpToReal(r), pressure_dev.interpToReal(r), emin, emax, lambda[0]); const Real sie = u / mass_density_dev.interpToReal(r); diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp index 156e069a0..1cc2a30c1 100644 --- a/src/radiation/cooling_function.cpp +++ b/src/radiation/cooling_function.cpp @@ -166,10 +166,10 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshBlockData *rc, const doub // Light Bulb with Liebendorfer model const bool is_monopole_cart = - (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleCart)); const bool is_monopole_sph = - (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph)); - + (typeid(PHOEBUS_GEOMETRY) == typeid(Geometry::MonopoleSph)); + auto &coords = pmb->coords; using Transformation_t = Geometry::SphericalToCartesian; auto const &geom_pkg = pmb->packages.Get("geometry"); @@ -189,14 +189,14 @@ TaskStatus CoolingFunctionCalculateFourForce(MeshBlockData *rc, const doub DEFAULT_LOOP_PATTERN, "CoolingFunctionCalculateFourForce", DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int k, const int j, const int i) { - const Real x1 = coords.Xc<1>(k, j, i); + const Real x1 = coords.Xc<1>(k, j, i); const Real x2 = coords.Xc<2>(k, j, i); const Real x3 = coords.Xc<3>(k, j, i); Real Cart[3]; Real s2c[3][3], c2s[3][3]; Real r; - if (is_monopole_sph){ + if (is_monopole_sph) { r = std::abs(x1); } else { r = std::sqrt(x1 * x1 + x2 * x2 + x3 * x3);