Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions inputs/radiation/crooked_pipe.in
Original file line number Diff line number Diff line change
Expand Up @@ -114,3 +114,4 @@ use_ddmc = false # use DDMC?
cutoff = 1.0e-6 # default is 1.0e-6
emit_temp_threshold = 1.2e4 # from ryan, 1.1e4 is background
#min_swarm_occupancy = 0.5
frequency_type = gray
1 change: 1 addition & 0 deletions inputs/radiation/rad_shock.in
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ temp_exp = 0.0 # temp exponent for opacity powerlaw
<radiation/imc>
num_particles = 100000 # particle resolution
use_ddmc = false # use DDMC?
frequency_type = gray

<problem>
rhol = 5.69 # density (left)
Expand Down
1 change: 1 addition & 0 deletions inputs/radiation/thermalization_imc.in
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ kappa_a = 2.0
<radiation/imc>
num_particles = 200000
use_ddmc = true
frequency_type = gray

# <radiation/moment>
# cfl = 0.3
Expand Down
2 changes: 2 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,8 @@ set (SRC_LIST
pgen/strat.hpp
pgen/thermalization.hpp

radiation/gas_opacity.cpp
radiation/gas_opacity.hpp
radiation/radiation.cpp
radiation/radiation.hpp
radiation/imc/imc_driver.cpp
Expand Down
23 changes: 18 additions & 5 deletions src/artemis.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -185,14 +185,27 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
// Operator split radiation
if (do_radiation) {
// Top-level radiation package
packages.Add(Radiation::Initialize(pin.get(), constants, do_imc));
packages.Add(Radiation::Initialize(pin.get(), units, constants, do_imc));
// Select between Jaybenne IMC or Moments
if (do_imc) {
auto eos_h = packages.Get("gas")->Param<EOS>("eos_h");
auto opacity_h = packages.Get("gas")->Param<MeanOpacity>("opacity_h");
auto scattering_h = packages.Get("gas")->Param<MeanScattering>("scattering_h");
packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h,
"radiation/imc"));
const auto frequency_type =
packages.Get("radiation")->Param<FrequencyType>("frequency_type");
if (frequency_type == FrequencyType::gray) {
auto opacity_h = packages.Get("radiation")->Param<MeanOpacity>("opacity_h");
auto scattering_h =
packages.Get("radiation")->Param<MeanScattering>("scattering_h");
packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h,
"radiation/imc"));
} else {
PARTHENON_REQUIRE(frequency_type == FrequencyType::multigroup,
"Invalid frequency_type!");
auto opacity_h = packages.Get("radiation")->Param<Opacity>("mg_opacity_h");
auto scattering_h =
packages.Get("radiation")->Param<Scattering>("mg_scattering_h");
packages.Add(jaybenne::Initialize(pin.get(), opacity_h, scattering_h, eos_h,
"radiation/imc"));
}
PARTHENON_REQUIRE(coords == Coordinates::cartesian,
"Jaybenne currently supports only Cartesian coordinates!");
} else if (do_moment) {
Expand Down
3 changes: 3 additions & 0 deletions src/artemis.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,9 @@ enum class ArtemisBC {
// Tensor indexing (currently used in radiation moments)
enum TensIdx { X11 = 0, X22 = 1, X33 = 2, X23 = 3, X13 = 4, X12 = 5 };

// Radiation solver frequency type enum
enum class FrequencyType { gray, multigroup };

// Floating point limits
template <typename T = Real>
KOKKOS_FORCEINLINE_FUNCTION constexpr auto Big() {
Expand Down
98 changes: 0 additions & 98 deletions src/gas/gas.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@
#include "utils/fluxes/fluid_fluxes.hpp"
#include "utils/history.hpp"
#include "utils/integrators/artemis_integrator.hpp"
#include "utils/opacity/opacity.hpp"
#include "utils/refinement/amr_criteria.hpp"
#include "utils/units.hpp"

Expand All @@ -43,7 +42,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
ArtemisUtils::Units &units,
ArtemisUtils::Constants &constants,
Packages_t &packages) {
using namespace singularity::photons;

auto gas = std::make_shared<StateDescriptor>("gas");
Params &params = gas->AllParams();
Expand Down Expand Up @@ -207,102 +205,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
}
params.Add("rsolver", riemann_solver);

// Opacity models
const Real time = units.GetTimeCodeToPhysical();
const Real mass = units.GetMassCodeToPhysical();
const Real length = units.GetLengthCodeToPhysical();
const Real temp = units.GetTemperatureCodeToPhysical();

// Absorption opacity model
std::string opacity_model_name =
pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant");

// Mean absorption opacity (either read from table or uses model
ArtemisUtils::MeanOpacity opacity;
if (opacity_model_name == "table") {
std::string table_filename =
pin->GetString("gas/opacity/absorption", "opacity_table");
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(table_filename), time, mass, length,
temp);
} else {
// Instantiate mean absorption opacity object (i.e., table)
const Real lRhoMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMin", -1.0);
const Real lRhoMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMax", 1.0);
const int NRho_a = pin->GetOrAddInteger("gas/opacity/absorption", "NRho", 2);
const Real lTMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMin", -1.0);
const Real lTMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMax", 1.0);
const int NT_a = pin->GetOrAddInteger("gas/opacity/absorption", "NT", 2);

if (opacity_model_name == "none") {
auto model = Gray(0.0);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
} else if (opacity_model_name == "constant") {
const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0);
auto model = Gray(kappa_a);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
} else if (opacity_model_name == "powerlaw") {
const Real coef_kappa_a =
pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0);
const Real rho_exp = pin->GetOrAddReal("gas/opacity/absorption", "rho_exp", 0.0);
const Real temp_exp = pin->GetOrAddReal("gas/opacity/absorption", "temp_exp", 0.0);
auto model = PowerLaw(coef_kappa_a, rho_exp, temp_exp);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
} else {
PARTHENON_FAIL("Opacity model not recognized!");
}
}
params.Add("opacity_h", opacity);
params.Add("opacity_d", opacity.GetOnDevice());

// Scattering opacity model
std::string scattering_model_name =
pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none");

// Instantiate mean scattering opacity object (i.e., table)
const Real lRhoMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMin", -1.0);
const Real lRhoMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMax", 1.0);
const int NRho_s = pin->GetOrAddInteger("gas/opacity/scattering", "NRho", 2);
const Real lTMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMin", -1.0);
const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0);
const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2);

ArtemisUtils::MeanScattering scattering;
if (scattering_model_name == "none") {
auto smodel = GrayS(0.0, 1.0);
scattering =
singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>(
singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s,
lTMin_s, lTMax_s, NT_s),
time, mass, length, temp);
} else if (scattering_model_name == "constant") {
const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0);
auto smodel = GrayS(kappa_s, 1.0);
scattering =
singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>(
singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s,
lTMin_s, lTMax_s, NT_s),
time, mass, length, temp);
} else {
PARTHENON_FAIL("Scattering model not recognized!");
}

params.Add("scattering_h", scattering);
params.Add("scattering_d", scattering.GetOnDevice());

// Dual energy switch
// When internal > de_switch * total we use the total
// The default turns off the switch
Expand Down
154 changes: 154 additions & 0 deletions src/radiation/gas_opacity.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
//========================================================================================
// (C) (or copyright) 2026. 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 "gas_opacity.hpp"
#include "utils/opacity/opacity.hpp"

namespace Gas {
void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units,
Params &params) {
using namespace singularity::photons;

// Opacity models
const Real time = units.GetTimeCodeToPhysical();
const Real mass = units.GetMassCodeToPhysical();
const Real length = units.GetLengthCodeToPhysical();
const Real temp = units.GetTemperatureCodeToPhysical();

// Get frequency type (it should already be set in radiation Initialization)
const auto frequency_type = params.Get<FrequencyType>("frequency_type");

// Absorption opacity model
std::string opacity_model_name =
pin->GetOrAddString("gas/opacity/absorption", "opacity_model", "constant");

// Mean absorption opacity (either read from table or uses model
ArtemisUtils::Opacity mg_opacity;
ArtemisUtils::MeanOpacity opacity;
if (opacity_model_name == "table") {
PARTHENON_REQUIRE(frequency_type == FrequencyType::gray,
"Only gray table opacity permitted, for now.");
std::string table_filename =
pin->GetString("gas/opacity/absorption", "opacity_table");
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(table_filename), time, mass, length,
temp);
} else {
// Instantiate mean absorption opacity object (i.e., table)
const Real lRhoMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMin", -1.0);
const Real lRhoMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lRhoMax", 1.0);
const int NRho_a = pin->GetOrAddInteger("gas/opacity/absorption", "NRho", 2);
const Real lTMin_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMin", -1.0);
const Real lTMax_a = pin->GetOrAddReal("gas/opacity/absorption", "lTMax", 1.0);
const int NT_a = pin->GetOrAddInteger("gas/opacity/absorption", "NT", 2);

if (opacity_model_name == "none") {
auto model = Gray(0.0);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
if (frequency_type == FrequencyType::multigroup) {
mg_opacity = singularity::photons::NonCGSUnits<singularity::photons::Gray>(
std::move(model), time, mass, length, temp);
}
} else if (opacity_model_name == "constant") {
const Real kappa_a = pin->GetOrAddReal("gas/opacity/absorption", "kappa_a", 0.0);
auto model = Gray(kappa_a);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
if (frequency_type == FrequencyType::multigroup) {
mg_opacity = singularity::photons::NonCGSUnits<singularity::photons::Gray>(
std::move(model), time, mass, length, temp);
}
} else if (opacity_model_name == "powerlaw") {
const Real coef_kappa_a =
pin->GetOrAddReal("gas/opacity/absorption", "coef_kappa_a", 0.0);
const Real rho_exp = pin->GetOrAddReal("gas/opacity/absorption", "rho_exp", 0.0);
const Real temp_exp = pin->GetOrAddReal("gas/opacity/absorption", "temp_exp", 0.0);
auto model = PowerLaw(coef_kappa_a, rho_exp, temp_exp);
opacity =
singularity::photons::MeanNonCGSUnits<singularity::photons::MeanOpacityBase>(
singularity::photons::MeanOpacityBase(model, lRhoMin_a, lRhoMax_a, NRho_a,
lTMin_a, lTMax_a, NT_a),
time, mass, length, temp);
if (frequency_type == FrequencyType::multigroup) {
mg_opacity = singularity::photons::NonCGSUnits<singularity::photons::PowerLaw>(
std::move(model), time, mass, length, temp);
}
} else {
PARTHENON_FAIL("Opacity model not recognized!");
}
}

params.Add("opacity_h", opacity);
params.Add("opacity_d", opacity.GetOnDevice());
if (frequency_type == FrequencyType::multigroup) {
params.Add("mg_opacity_h", mg_opacity);
params.Add("mg_opacity_d", mg_opacity.GetOnDevice());
}

// Scattering opacity model
std::string scattering_model_name =
pin->GetOrAddString("gas/opacity/scattering", "scattering_model", "none");

// Instantiate mean scattering opacity object (i.e., table)
const Real lRhoMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMin", -1.0);
const Real lRhoMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lRhoMax", 1.0);
const int NRho_s = pin->GetOrAddInteger("gas/opacity/scattering", "NRho", 2);
const Real lTMin_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMin", -1.0);
const Real lTMax_s = pin->GetOrAddReal("gas/opacity/scattering", "lTMax", 1.0);
const int NT_s = pin->GetOrAddInteger("gas/opacity/scattering", "NT", 2);

ArtemisUtils::Scattering mg_scattering;
ArtemisUtils::MeanScattering scattering;
if (scattering_model_name == "none") {
auto smodel = GrayS(0.0, 1.0);
scattering =
singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>(
singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s,
lTMin_s, lTMax_s, NT_s),
time, mass, length, temp);
if (frequency_type == FrequencyType::multigroup) {
mg_scattering = singularity::photons::NonCGSUnitsS<singularity::photons::GrayS>(
std::move(smodel), time, mass, length, temp);
}
} else if (scattering_model_name == "constant") {
const Real kappa_s = pin->GetOrAddReal("gas/opacity/scattering", "kappa_s", 0.0);
auto smodel = GrayS(kappa_s, 1.0);
scattering =
singularity::photons::MeanNonCGSUnitsS<singularity::photons::MeanSOpacityCGS>(
singularity::photons::MeanSOpacityCGS(smodel, lRhoMin_s, lRhoMax_s, NRho_s,
lTMin_s, lTMax_s, NT_s),
time, mass, length, temp);
if (frequency_type == FrequencyType::multigroup) {
mg_scattering = singularity::photons::NonCGSUnitsS<singularity::photons::GrayS>(
std::move(smodel), time, mass, length, temp);
}
} else {
PARTHENON_FAIL("Scattering model not recognized!");
}

params.Add("scattering_h", scattering);
params.Add("scattering_d", scattering.GetOnDevice());
if (frequency_type == FrequencyType::multigroup) {
params.Add("mg_scattering_h", mg_scattering);
params.Add("mg_scattering_d", mg_scattering.GetOnDevice());
}
}

} // namespace Gas
28 changes: 28 additions & 0 deletions src/radiation/gas_opacity.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
//========================================================================================
// (C) (or copyright) 2026. 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.
//========================================================================================
#ifndef RADIATION_GAS_OPACITY_HPP_
#define RADIATION_GAS_OPACITY_HPP_

// Artemis includes
#include "utils/units.hpp"

// Parthenon includes
#include <parthenon/driver.hpp>
#include <parthenon/package.hpp>

namespace Gas {
void InitGasOpacity(ParameterInput *pin, const ArtemisUtils::Units &units,
Params &params);
} // namespace Gas

#endif // RADIATION_GAS_OPACITY_HPP_
Loading
Loading