Skip to content

Commit

Permalink
debugging
Browse files Browse the repository at this point in the history
  • Loading branch information
carlnotsagan committed Jul 5, 2023
1 parent d2172ae commit 3906a33
Show file tree
Hide file tree
Showing 2 changed files with 42 additions and 69 deletions.
34 changes: 12 additions & 22 deletions inputs/standing_accretion_shock.pin
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -79,10 +79,15 @@ field = c.c.bulk.rho
method = derivative_order_1
max_level = 1

#<eos>
#type = IdealGas
#Gamma = 1.33
#Cv = 3.0

<eos>
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

<physics>
hydro = true
Expand All @@ -92,7 +97,7 @@ rad = false

<fluid>
xorder = 2
cfl = 0.5
cfl = 0.6
riemann = llf
recon = weno5
c2p_method = robust
Expand All @@ -101,24 +106,9 @@ c2p_tol = 1.e-8
Ye = true
mhd = false

<fixup>
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


<geometry>
a = 0

<coordinates>
derefine_poles = false

<units>
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
Expand Down
77 changes: 30 additions & 47 deletions src/pgen/standing_accretion_shock.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
// Parthenon
#include <globals.hpp>

#include "geometry/boyer_lindquist.hpp"
#include "pgen/pgen.hpp"
#include "phoebus_utils/root_find.hpp"
#include "phoebus_utils/unit_conversions.hpp"
Expand Down Expand Up @@ -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(
Expand All @@ -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);
Expand All @@ -98,48 +97,37 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
auto eos = pmb->packages.Get("eos")->Param<Microphysics::EOS::EOS>("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<Real>("Cv");
auto epsmin = pmb->packages.Get("eos")->Param<Real>("sie_min");
auto epsmax = pmb->packages.Get("eos")->Param<Real>("sie_max");
//auto Cv = pmb->packages.Get("eos")->Param<Real>("Cv");
auto Tmin = pmb->packages.Get("eos")->Param<Real>("T_min");
auto Tmax = pmb->packages.Get("eos")->Param<Real>("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<bool>("derefine_poles");
Real h = gpkg->Param<Real>("h");
Real xt = gpkg->Param<Real>("xt");
Real alpha = gpkg->Param<Real>("alpha");
Real x0 = gpkg->Param<Real>("x0");
Real smooth = gpkg->Param<Real>("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) {
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 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;
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];
Real eos_lambda[2] = {0.};
if (iye > 0) {
v(iye, k, j, i) = 0.5;
eos_lambda[0] = v(iye, k, j, i);
Expand All @@ -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.;
Expand All @@ -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];
Expand All @@ -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);
Expand All @@ -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
Expand All @@ -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;
}

Expand Down

0 comments on commit 3906a33

Please sign in to comment.