diff --git a/inputs/standing_accretion_shock.pin b/inputs/standing_accretion_shock.pin new file mode 100644 index 00000000..b25e6c3e --- /dev/null +++ b/inputs/standing_accretion_shock.pin @@ -0,0 +1,145 @@ +# © 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 = reflect +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, & + temperature,& + p.ye, & + g.c.coord, & + g.n.coord, & + flux_divergence, & + src_terms + + +file_type = hdf5 # Tabular data dump +dt = 100 # time increment between outputs + + +nlim = -1 # cycle limit +tlim = 10000 # time limit +integrator = rk2 # time integration algorithm +ncycle_out = 100 # interval for stdout summary info +dt_init_fact = 1.e-8 + + +nghost = 4 +#refinement = adaptive +#numlevel = 3 + +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 +ox1_bc = user # Outer-X1 boundary condition flag + +nx2 = 1 # Number of zones in X2-direction +x2min = 0 # minimum value of X2 +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 + +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 = 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 + + +enable_floors = false +enable_mhd_floors = false +enable_ceilings = false +enable_flux_fixup = false +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 +he = false +3t = false +rad = true + + +xorder = 2 +cfl = 0.6 +riemann = llf +recon = weno5 +c2p_method = robust +c2p_max_iter = 1000 +c2p_tol = 1.e-8 +Ye = true +mhd = false + + +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 = 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_10k.txt + + +method=cooling_function +do_liebendorfer=false +do_lightbulb=true +lum=2.5 +do_gain_reducer=false +always_gain=false diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5fa94e05..6737e408 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 2df82523..048edd0c 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 \ diff --git a/src/pgen/progenitor.cpp b/src/pgen/progenitor.cpp index 0c6cda5d..b6409e00 100644 --- a/src/pgen/progenitor.cpp +++ b/src/pgen/progenitor.cpp @@ -124,7 +124,6 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { 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), pressure_dev.interpToReal(r), emin, emax, lambda[0]); diff --git a/src/pgen/standing_accretion_shock.cpp b/src/pgen/standing_accretion_shock.cpp new file mode 100644 index 00000000..7dcff2ab --- /dev/null +++ b/src/pgen/standing_accretion_shock.cpp @@ -0,0 +1,259 @@ +// 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. +// 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. + +// 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; + +constexpr int NSASW = 8; +constexpr int R = 0; +constexpr int RHO = 1; +constexpr int VR = 2; +constexpr int EPS = 3; + +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); + +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!"); + + auto &rc = pmb->meshblock_data.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; + + const std::string eos_type = pin->GetString("eos", "type"); + PARTHENON_REQUIRE_THROWS( + eos_type == "IdealGas" || eos_type == "StellarCollapse", + "Standing Accretion Shock setup only works with Ideal Gas or Stellar Collapse EOS"); + + 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"); + const Real a = pin->GetReal("geometry", "a"); + auto geom = Geometry::GetCoordinateSystem(rc.get()); + + // info about model file + std::string model_filename = + pin->GetOrAddString("standing_accretion_shock", "model_filename", "file.txt"); + + 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, NSASW, state_raw_h); + + Kokkos::deep_copy(state_raw_d, state_raw_h); + + const Real rhomin = 1.e-16; + const Real epsmin = 1.e-16; + + 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); + + Radius raw_radius(rin, rout, num_zones); + + 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) { + 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 = std::abs(x1); + Real eos_lambda[2] = {0.}; + + if (iye > 0) { + v(iye, k, j, i) = 0.5; + eos_lambda[0] = v(iye, k, j, i); + } + + 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)); + + v(irho, k, j, i) = rho; + v(ieng, k, j, i) = rho * eps; + v(iprs, k, j, i) = eos.PressureFromDensityInternalEnergy(rho, eps, eos_lambda); + v(igm1, k, j, i) = + eos.BulkModulusFromDensityInternalEnergy(rho, eps, eos_lambda) / + v(iprs, k, j, i); + + // 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); + + for (int d = 0; d < 3; d++) { + v(ivlo + d, k, j, i) = ucon[d + 1]; + } + }); + 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) { + + 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 +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 diff --git a/src/radiation/cooling_function.cpp b/src/radiation/cooling_function.cpp index bd47ddbe..1cc2a30c 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) { @@ -181,7 +189,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