Skip to content

Precompute gray opacity #78

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
adamdempsey90 merged 11 commits into from
Jun 26, 2025
Merged

Precompute gray opacity #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
5809482
Precompute gray opacity per cell.
RyanWollaeger Jun 19, 2025
242c540
Move opacity precomputation to dedicated task (collection).
RyanWollaeger Jun 19, 2025
a4cec01
@pdmullen pass at A--J API change
pdmullen Jun 20, 2025
b2e2ec6
Opt in to compiler timings
pdmullen Jun 20, 2025
d64fa64
Consistency in ordering
pdmullen Jun 20, 2025
14399b4
Bump jaybenne
pdmullen Jun 20, 2025
9d13373
Merge remote-tracking branch 'origin/develop' into precomp_gray_opac
pdmullen Jun 20, 2025
d73a7f4
Extract params from radiation pkg
pdmullen Jun 20, 2025
5e71ed0
Bump jaybenne
pdmullen Jun 25, 2025
88388f4
Bump jaybenne again
pdmullen Jun 25, 2025
b8e8d64
Bump jaybenne
pdmullen Jun 25, 2025
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
4 changes: 3 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ option(ARTEMIS_ENABLE_HIP "Enable hip for artemis and all dependencies" OFF)
option(ARTEMIS_ENABLE_HDF5 "Enable HDF5 for artemis and all dependencies" ON)
option(ARTEMIS_ENABLE_MPI "Enable MPI for artemis and all dependencies" ON)
option(ARTEMIS_ENABLE_OPENMP "Enable OpenMP for artemis and parthenon" OFF)
option(ARTEMIS_ENABLE_COMPILE_TIMING "Enable timing of compilation of artemis" ON)
option(ARTEMIS_ENABLE_COMPILE_TIMING "Enable timing of compilation of artemis" OFF)
option(ARTEMIS_ENABLE_ASAN "Enable AddressSanitizer to detect memory errors" OFF)
option(ARTEMIS_PORTABLE_RESTART "Enable portable data types for REBOUND restarts" ON)

Expand Down Expand Up @@ -228,6 +228,8 @@ set(JAYBENNE_HOST_VARIABLE_HEADER "src/artemis.hpp")
set(JAYBENNE_HOST_DENSITY_VARIABLE "gas::prim::density")
set(JAYBENNE_HOST_SPECIFIC_INTERNAL_ENERGY_VARIABLE "gas::prim::sie")
set(JAYBENNE_HOST_UPDATE_ENERGY_VARIABLE "gas::cons::internal_energy")
set(JAYBENNE_HOST_ABSORPTION_OPACITY_VARIABLE "rad::opac::absorption")
set(JAYBENNE_HOST_SCATTERING_OPACITY_VARIABLE "rad::opac::scattering")

# Add jaybenne
message(STATUS "Adding jaybenne and dependencies")
Expand Down
2 changes: 1 addition & 1 deletion external/jaybenne
Open in desktop
Submodule jaybenne updated 10 files
+1 −1 external/singularity-eos
+6 −0 src/CMakeLists.txt
+3 −0 src/jaybenne/jaybenne.cpp
+2 −0 src/jaybenne/jaybenne_config.hpp.in
+2 −0 src/jaybenne/jaybenne_variables.hpp
+2 −1 src/jaybenne/sourcing.cpp
+9 −14 src/jaybenne/transport.cpp
+5 −11 src/jaybenne/transport_ddmc.cpp
+66 −7 src/mcblock/mcblock.cpp
+2 −0 src/mcblock/mcblock_variables.hpp
2 changes: 2 additions & 0 deletions src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ set (SRC_LIST
pgen/strat.hpp
pgen/thermalization.hpp

radiation/radiation.cpp
radiation/radiation.hpp
radiation/imc/imc_driver.cpp
radiation/imc/imc.hpp
radiation/moments/moments.cpp
Expand Down
7 changes: 5 additions & 2 deletions src/artemis.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "gravity/gravity.hpp"
#include "nbody/nbody.hpp"
#include "radiation/moments/moments.hpp"
#include "radiation/radiation.hpp"
#include "rotating_frame/rotating_frame.hpp"
#include "utils/artemis_utils.hpp"
#include "utils/history.hpp"
Expand Down Expand Up @@ -127,7 +128,9 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
if (do_cooling) packages.Add(Gas::Cooling::Initialize(pin.get()));
if (do_drag) packages.Add(Drag::Initialize(pin.get()));
if (do_radiation) {
// swap between native artemis radiation and jaybenne imc
// Top-level radiation package
packages.Add(Radiation::Initialize(pin.get(), 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");
Expand All @@ -137,7 +140,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
PARTHENON_REQUIRE(coords == Coordinates::cartesian,
"Jaybenne currently supports only Cartesian coordinates!");
} else if (do_moment) {
packages.Add(Radiation::Initialize(pin.get(), constants));
packages.Add(Moments::Initialize(pin.get(), constants));
} else {
PARTHENON_FAIL("Unknown radiation model!");
}
Expand Down
4 changes: 4 additions & 0 deletions src/artemis.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,10 @@ ARTEMIS_VARIABLE(rad.prim, energy);
ARTEMIS_VARIABLE(rad.prim, pressure);
ARTEMIS_VARIABLE(rad.prim, flux);
} // namespace prim
namespace opac {
ARTEMIS_VARIABLE(rad.opac, absorption);
ARTEMIS_VARIABLE(rad.opac, scattering);
} // namespace opac
} // namespace rad

#undef ARTEMIS_VARIABLE
Expand Down
2 changes: 1 addition & 1 deletion src/artemis_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ TaskListStatus ArtemisDriver<GEOM>::Step() {
if (status != TaskListStatus::complete) return status;

// Operator split, moments subcyling (M1 or P1)
if (do_moment) status = Radiation::MomentsDriver<GEOM>(pmesh, tm, rad_integrator.get());
if (do_moment) status = Moments::MomentsDriver<GEOM>(pmesh, tm, rad_integrator.get());
if (status != TaskListStatus::complete) return status;

// Compute new dt, (de)refine, and handle sparse (if enabled)
Expand Down
8 changes: 3 additions & 5 deletions src/derived/fill_derived.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@ namespace ArtemisDerived {
template <Coordinates GEOM>
TaskStatus SetAuxillaryFields(MeshData<Real> *md) {
using parthenon::MakePackDescriptor;
using TE = parthenon::TopologicalElement;
auto pm = md->GetParentPointer();
auto &resolved_pkgs = pm->resolved_packages;

Expand Down Expand Up @@ -84,7 +83,6 @@ TaskStatus SetAuxillaryFields(MeshData<Real> *md) {
template <Coordinates GEOM>
void ConsToPrim(MeshData<Real> *md) {
using parthenon::MakePackDescriptor;
using TE = parthenon::TopologicalElement;
auto pm = md->GetParentPointer();
auto &resolved_pkgs = pm->resolved_packages;

Expand Down Expand Up @@ -193,7 +191,7 @@ void ConsToPrim(MeshData<Real> *md) {
const Real hfx1 = vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) / conv[0];
const Real hfx2 = vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) / conv[1];
const Real hfx3 = vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) / conv[2];
const auto fx = Radiation::NormalizeFlux(hfx1, hfx2, hfx3);
const auto fx = Moments::NormalizeFlux(hfx1, hfx2, hfx3);
vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0];
vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1];
vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2];
Expand All @@ -211,7 +209,6 @@ void ConsToPrim(MeshData<Real> *md) {
template <typename T, Coordinates GEOM>
void PrimToCons(T *md) {
using parthenon::MakePackDescriptor;
using TE = parthenon::TopologicalElement;
auto pm = md->GetParentPointer();
auto &resolved_pkgs = pm->resolved_packages;

Expand All @@ -220,6 +217,7 @@ void PrimToCons(T *md) {
const bool do_gas = artemis_pkg->template Param<bool>("do_gas");
const bool do_dust = artemis_pkg->template Param<bool>("do_dust");
const bool do_rad = artemis_pkg->template Param<bool>("do_moment");
const bool do_imc = artemis_pkg->template Param<bool>("do_imc");

// Extract gas parameters
Real dflr_gas = Null<Real>();
Expand Down Expand Up @@ -344,7 +342,7 @@ void PrimToCons(T *md) {
const Real fx1 = vmesh(b, rad::prim::flux(VI(n, 0)), k, j, i);
const Real fx2 = vmesh(b, rad::prim::flux(VI(n, 1)), k, j, i);
const Real fx3 = vmesh(b, rad::prim::flux(VI(n, 2)), k, j, i);
const auto fx = Radiation::NormalizeFlux(fx1, fx2, fx3);
const auto fx = Moments::NormalizeFlux(fx1, fx2, fx3);
vmesh(b, rad::cons::flux(VI(n, 0)), k, j, i) = fx[0] * conv[0];
vmesh(b, rad::cons::flux(VI(n, 1)), k, j, i) = fx[1] * conv[1];
vmesh(b, rad::cons::flux(VI(n, 2)), k, j, i) = fx[2] * conv[2];
Expand Down
5 changes: 4 additions & 1 deletion src/radiation/imc/imc_driver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include "artemis.hpp"
#include "derived/fill_derived.hpp"
#include "radiation/imc/imc.hpp"
#include "radiation/radiation.hpp"

// Jaybenne includes
#include "jaybenne.hpp"
Expand All @@ -28,7 +29,9 @@ namespace IMC {
//! \brief Executes thermal IMC transport (Jaybenne) and syncs updated fields
template <Coordinates GEOM>
TaskListStatus JaybenneIMC(Mesh *pmesh, const Real time, const Real dt) {
auto status = jaybenne::RadiationStep(pmesh, time, dt).Execute();
auto status = Radiation::UpdateRadiationFields(pmesh).Execute();
if (status != TaskListStatus::complete) return status;
status = jaybenne::RadiationStep(pmesh, time, dt).Execute();
if (status != TaskListStatus::complete) return status;
status = ArtemisDerived::SyncFields<GEOM>(pmesh, time, dt).Execute();
return status;
Expand Down
48 changes: 24 additions & 24 deletions src/radiation/moments/matter_coupling.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,10 @@ using ArtemisUtils::MeanOpacity;
using ArtemisUtils::MeanScattering;
using ArtemisUtils::VI;

namespace Radiation {
namespace Moments {

//----------------------------------------------------------------------------------------
//! \fn TaskStatus Radiation::MatterCouplingSimpleImpl
//! \fn TaskStatus Moments::MatterCouplingSimpleImpl
//! \brief Implementation for simple radiation-matter coupling source
template <Coordinates GEOM, Closure CLOSURE>
TaskStatus MatterCouplingSimpleImpl(MeshData<Real> *u0, const Real dt) {
Expand All @@ -47,15 +47,15 @@ TaskStatus MatterCouplingSimpleImpl(MeshData<Real> *u0, const Real dt) {
auto dflr = gas_pkg->template Param<Real>("dfloor");
auto de_switch = gas_pkg->template Param<Real>("de_switch");

// Extract radiation package and params
auto &rad_pkg = pm->packages.Get("moments");
const auto chat = rad_pkg->template Param<Real>("chat");
const auto c = rad_pkg->template Param<Real>("c");
const auto arad = rad_pkg->template Param<Real>("arad");
const auto outer_max = rad_pkg->template Param<int>("outer_iteration_max");
const auto inner_max = rad_pkg->template Param<int>("inner_iteration_max");
const auto outer_tol = rad_pkg->template Param<Real>("outer_iteration_tol");
const auto inner_tol = rad_pkg->template Param<Real>("inner_iteration_tol");
// Extract radiation and moments package and params
auto &moments_pkg = pm->packages.Get("moments");
const auto chat = moments_pkg->template Param<Real>("chat");
const auto c = moments_pkg->template Param<Real>("c");
const auto arad = moments_pkg->template Param<Real>("arad");
const auto outer_max = moments_pkg->template Param<int>("outer_iteration_max");
const auto inner_max = moments_pkg->template Param<int>("inner_iteration_max");
const auto outer_tol = moments_pkg->template Param<Real>("outer_iteration_tol");
const auto inner_tol = moments_pkg->template Param<Real>("inner_iteration_tol");

// Extract rotating frame quantities
Real om0 = 0.0;
Expand All @@ -80,7 +80,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData<Real> *u0, const Real dt) {
// Prepare scratch pad memory
// const int ncells1 = ib.e - ib.s + 1 + 2 * parthenon::Globals::nghost;
// int scr_size = ScratchPad1D<Real>::shmem_size(ncells1) * 12;
// const int scr_level = rad_pkg->template Param<int>("scr_level");
// const int scr_level = moments_pkg->template Param<int>("scr_level");
parthenon::par_for(
DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1,
kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
Expand Down Expand Up @@ -182,7 +182,7 @@ TaskStatus MatterCouplingSimpleImpl(MeshData<Real> *u0, const Real dt) {
}

//----------------------------------------------------------------------------------------
//! \fn TaskStatus Radiation::MatterCouplingSimpleImpl
//! \fn TaskStatus Moments::MatterCouplingSimpleImpl
//! \brief Implementation for "full" radiation-matter coupling source
template <Coordinates GEOM, Closure CLOSURE>
TaskStatus MatterCouplingFullSingleImpl(MeshData<Real> *u0, const Real dt) {
Expand All @@ -201,14 +201,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData<Real> *u0, const Real dt) {
auto de_switch = gas_pkg->template Param<Real>("de_switch");

// Extract radiation package and params
auto &rad_pkg = pm->packages.Get("moments");
const auto chat = rad_pkg->template Param<Real>("chat");
const auto c = rad_pkg->template Param<Real>("c");
const auto arad = rad_pkg->template Param<Real>("arad");
const auto outer_max = rad_pkg->template Param<int>("outer_iteration_max");
const auto inner_max = rad_pkg->template Param<int>("inner_iteration_max");
const auto outer_tol = rad_pkg->template Param<Real>("outer_iteration_tol");
const auto inner_tol = rad_pkg->template Param<Real>("inner_iteration_tol");
auto &moments_pkg = pm->packages.Get("moments");
const auto chat = moments_pkg->template Param<Real>("chat");
const auto c = moments_pkg->template Param<Real>("c");
const auto arad = moments_pkg->template Param<Real>("arad");
const auto outer_max = moments_pkg->template Param<int>("outer_iteration_max");
const auto inner_max = moments_pkg->template Param<int>("inner_iteration_max");
const auto outer_tol = moments_pkg->template Param<Real>("outer_iteration_tol");
const auto inner_tol = moments_pkg->template Param<Real>("inner_iteration_tol");

// Extract rotating frame quantities
Real om0 = 0.0;
Expand All @@ -234,7 +234,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData<Real> *u0, const Real dt) {
// Prepare scratch pad memory
// const int ncells1 = ib.e - ib.s + 1 + 2 * parthenon::Globals::nghost;
// int scr_size = ScratchPad1D<Real>::shmem_size(ncells1) * 12;
// const int scr_level = rad_pkg->template Param<int>("scr_level");
// const int scr_level = moments_pkg->template Param<int>("scr_level");
parthenon::par_for(
DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1,
kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
Expand Down Expand Up @@ -421,6 +421,6 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData<Real> *u0, const Real dt) {
return TaskStatus::complete;
}

} // namespace Radiation
} // namespace Moments

#endif // RADIATION_MOMENTS_MATTER_COUPLING_HPP_
#endif // RADIATION_MOMENTS_MATTER_COUPLING_HPP_
Loading