Skip to content
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.

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

Gray vs multigroup #5

Open
wants to merge 17 commits into
base: develop
Choose a base branch
from
Open
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
5 changes: 2 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@ Currently supported computers/partitions are:

## Darwin

power9-rhel7
volta-x86
skylake-gold

Expand Down Expand Up @@ -92,5 +91,5 @@ of `jaybenne`.

# Run driver executable

cd build/src
mpiexec -n 1 ./mcblock -i ../../prob/jbinput.stepdiff
cd build
mpiexec -n 1 ./mcblock -i ../inputs/stepdiff.in
1 change: 1 addition & 0 deletions inputs/inf.in
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = constant
opacity_constant_value = 1.0 # cm^2/g
scattering_model = constant
Expand Down
1 change: 1 addition & 0 deletions inputs/inf_stiff.in
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ source_strategy = uniform
seed = 349856

<mcblock>
frequency_type = gray
opacity_model = constant
opacity_constant_value = 1000.0 # cm^2/g
scattering_model = constant
Expand Down
1 change: 1 addition & 0 deletions inputs/stepdiff.in
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = none
scattering_model = constant
scattering_constant_value = 1.0e3
Expand Down
1 change: 1 addition & 0 deletions inputs/stepdiff_ddmc.in
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = none
scattering_model = constant
scattering_constant_value = 1.0e3
Expand Down
1 change: 1 addition & 0 deletions inputs/stepdiff_smr.in
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = none
scattering_model = constant
scattering_constant_value = 1.0e3
Expand Down
1 change: 1 addition & 0 deletions inputs/stepdiff_smr_ddmc.in
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = none
scattering_model = constant
scattering_constant_value = 1.0e3
Expand Down
1 change: 1 addition & 0 deletions inputs/stepdiff_smr_hybrid.in
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ source_strategy = uniform
seed = 349857

<mcblock>
frequency_type = gray
opacity_model = none
scattering_model = constant
scattering_constant_value = 1.0e3
Expand Down
6 changes: 6 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,12 @@ if(JAYBENNE_STANDALONE_MODE)
set(JAYBENNE_SCATTERING_OPACITY_TYPE
"mcblock::Scattering"
CACHE STRING "Namespace for scattering opacity Singularity type" FORCE)
set(JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE
"mcblock::MeanOpacity"
CACHE STRING "Namespace for absorption mean opacity Singularity type" FORCE)
set(JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE
"mcblock::MeanScattering"
CACHE STRING "Namespace for scattering mean opacity Singularity type" FORCE)
set(JAYBENNE_HOST_VARIABLE_HEADER
"../mcblock/mcblock_variables.hpp"
CACHE STRING "Path to host code variables header file" FORCE)
Expand Down
244 changes: 197 additions & 47 deletions src/jaybenne/jaybenne.cpp

Large diffs are not rendered by default.

9 changes: 8 additions & 1 deletion src/jaybenne/jaybenne.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,17 +50,24 @@ namespace jaybenne {
std::shared_ptr<parthenon::StateDescriptor> Initialize(parthenon::ParameterInput *pin,
Opacity &opacity,
Scattering &scattering, EOS &eos);
std::shared_ptr<parthenon::StateDescriptor> Initialize(parthenon::ParameterInput *pin,
MeanOpacity &mopacity,
MeanScattering &mscattering,
EOS &eos);

// Model enums
enum class SourceStrategy { uniform, energy };
enum class SourceType { thermal, emission };
enum class FrequencyType { gray, multigroup };

// Tasks
template <FrequencyType FT>
TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real dt);
template <FrequencyType FT>
TaskStatus TransportPhotons_DDMC(MeshData<Real> *md, const Real t_start, const Real dt);
TaskStatus SampleDDMCBlockFace(MeshData<Real> *md);
TaskStatus CheckCompletion(MeshData<Real> *md, const Real t_end);
template <typename T, SourceType ST>
template <typename T, SourceType ST, FrequencyType FT>
TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt);
TaskStatus DefragParticles(MeshBlock *pmb);
TaskStatus UpdateDerivedTransportFields(MeshData<Real> *md, const Real dt);
Expand Down
2 changes: 2 additions & 0 deletions src/jaybenne/jaybenne_config.hpp.in
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@
using @JAYBENNE_EOS_TYPE@;
using @JAYBENNE_ABSORPTION_OPACITY_TYPE@;
using @JAYBENNE_SCATTERING_OPACITY_TYPE@;
using @JAYBENNE_ABSORPTION_MEAN_OPACITY_TYPE@;
using @JAYBENNE_SCATTERING_MEAN_OPACITY_TYPE@;

#include "@JAYBENNE_HOST_VARIABLE_HEADER@"

Expand Down
1 change: 1 addition & 0 deletions src/jaybenne/jaybenne_variables.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ JAYBENNE_FIELD_VARIABLE(field.jaybenne, fleck_factor);
JAYBENNE_FIELD_VARIABLE(field.jaybenne, ddmc_face_prob);
JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_ew_per_cell);
JAYBENNE_FIELD_VARIABLE(field.jaybenne, source_num_per_cell);
JAYBENNE_FIELD_VARIABLE(field.jaybenne, emission_cdf);
JAYBENNE_FIELD_VARIABLE(field.jaybenne, energy_delta);
namespace host {
typedef HOST_DENSITY density;
Expand Down
104 changes: 93 additions & 11 deletions src/jaybenne/sourcing.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ namespace jaybenne {
//! \fn TaskStatus SourcePhotons
//! \brief Create photons, either during initialization or during a timestep.
//! TODO(BRR) modify interface so we don't need t_start, dt for initialization
template <typename T, SourceType ST>
template <typename T, SourceType ST, FrequencyType FT>
TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
namespace fj = field::jaybenne;
namespace fjh = field::jaybenne::host;
Expand All @@ -31,8 +31,21 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
auto pm = md->GetParentPointer();
auto &resolved_pkgs = pm->resolved_packages;
auto &jb_pkg = pm->packages.Get("jaybenne");
const Real h = jb_pkg->template Param<Real>("planck_constant");
auto &eos = jb_pkg->template Param<EOS>("eos_d");
auto &opacity = jb_pkg->template Param<Opacity>("opacity_d");
int n_nubins = -1;
Real numin = -1.;
Real numax = -1.;
MeanOpacity mopacity;
Opacity opacity;
if constexpr (FT == FrequencyType::gray) {
mopacity = jb_pkg->template Param<MeanOpacity>("mopacity_d");
} else if constexpr (FT == FrequencyType::multigroup) {
opacity = jb_pkg->template Param<Opacity>("opacity_d");
n_nubins = jb_pkg->template Param<int>("n_nubins");
numin = jb_pkg->template Param<Real>("numin");
numax = jb_pkg->template Param<Real>("numax");
}
RyanWollaeger marked this conversation as resolved.
Show resolved Hide resolved
auto &do_emission = jb_pkg->template Param<bool>("do_emission");
auto &source_strategy = jb_pkg->template Param<SourceStrategy>("source_strategy");
PARTHENON_REQUIRE(source_strategy != SourceStrategy::energy,
Expand All @@ -51,7 +64,8 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
// Create pack
static auto desc =
MakePackDescriptor<fjh::density, fjh::sie, fj::fleck_factor, fj::source_ew_per_cell,
fj::source_num_per_cell, fj::energy_delta>(resolved_pkgs.get());
fj::source_num_per_cell, fj::emission_cdf, fj::energy_delta>(
resolved_pkgs.get());
auto vmesh = desc.GetPack(md);

// Indexing and dimensionality
Expand Down Expand Up @@ -87,12 +101,47 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
[[maybe_unused]] const Real &sbd = sb;
[[maybe_unused]] const Real &vvd = vv;
[[maybe_unused]] const Real &dtd = dt;
[[maybe_unused]] auto &opac = opacity;
[[maybe_unused]] auto opac = opacity;
[[maybe_unused]] auto mopac = mopacity;
[[maybe_unused]] const auto numind = numin;
[[maybe_unused]] const auto numaxd = numax;
[[maybe_unused]] const auto n_nubinsd = n_nubins;
Real erad = 0.0;
if constexpr (ST == SourceType::thermal) {
erad = (4.0 * sbd / vvd) * std::pow(temp, 4.0) * dv;
} else if constexpr (ST == SourceType::emission) {
const Real emis = opac.Emissivity(rho, temp);
Real emis = -1.;
if constexpr (FT == FrequencyType::gray) {
emis = mopac.Emissivity(rho, temp);
} else if constexpr (FT == FrequencyType::multigroup) {
// Construct emission CDF
const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd;
// this is the mid-point of the 1st group in log-space:
// nu=exp(log(numin) + 0.5*dlnu)
Real nu = numind * std::exp(0.5 * dlnu);
Real dnu = nu * dlnu;
RyanWollaeger marked this conversation as resolved.
Show resolved Hide resolved
vmesh(b, fj::emission_cdf(0), k, j, i) =
opac.EmissivityPerNu(rho, temp, numind * std::exp(0.5 * dlnu)) *
dnu;
for (int n = 1; n < n_nubinsd; n++) {
// this is the mid-point of group n in log-space:
// nu=exp(log(numin)+(n+0.5)*dlnu)
nu = numind * std::exp((n + 0.5) * dlnu);
dnu = dlnu * nu;
RyanWollaeger marked this conversation as resolved.
Show resolved Hide resolved
vmesh(b, fj::emission_cdf(n), k, j, i) =
opac.EmissivityPerNu(rho, temp, nu) * dnu +
vmesh(b, fj::emission_cdf(n - 1), k, j, i);
}
emis = 0.;
for (int n = 0; n < n_nubinsd; n++) {
const Real dnu = dlnu * numind * std::exp((n + 0.5) * dlnu);
// Get total emissivity
emis += vmesh(b, fj::emission_cdf(n), k, j, i);
brryan marked this conversation as resolved.
Show resolved Hide resolved
// Normalize emission CDF
vmesh(b, fj::emission_cdf(n), k, j, i) /=
vmesh(b, fj::emission_cdf(n_nubinsd - 1), k, j, i);
}
}
erad = vmesh(b, fj::fleck_factor(), k, j, i) * emis * dv * dtd;
}
// Set source_num_per_cell
Expand Down Expand Up @@ -156,6 +205,10 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
auto rng_gen = rng_pool.get_state();
[[maybe_unused]] const Real &dtd = dt;
[[maybe_unused]] const Real &t_startd = t_start;
[[maybe_unused]] const Real hd = h;
[[maybe_unused]] const Real numaxd = numax;
[[maybe_unused]] const Real numind = numin;
[[maybe_unused]] const int n_nubinsd = n_nubins;

// Starting index and length of particles in this cell
const int &pstart_idx = prefix_sum(b, cell_idx_1d);
Expand All @@ -170,6 +223,9 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
ppack_i(b, ph::ijk(1), n) = j;
ppack_i(b, ph::ijk(2), n) = k;

// Set energy weight
ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i);

// Sample position uniformly in space over cell
// TODO(BRR) only valid for Cartesian
ppack_r(b, swarm_position::x(), n) = xi + dx_i * (rng_gen.drand() - 0.5);
Expand All @@ -189,8 +245,21 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
const Real &rho = vmesh(b, fjh::density(), k, j, i);
const Real &sie = vmesh(b, fjh::sie(), k, j, i);
const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp);
ppack_r(b, ph::weight(), n) = vmesh(b, fj::source_ew_per_cell(), k, j, i);
if constexpr (FT == FrequencyType::gray) {
ppack_r(b, ph::energy(), n) = sample_Planck_energy(rng_gen, sb, temp);
} else if constexpr (FT == FrequencyType::multigroup) {
// Sample energy from CDF
const Real rand = rng_gen.drand();
int n;
for (n = 0; n < n_nubinsd; n++) {
if (vmesh(b, fj::emission_cdf(n), k, j, i) >= rand) {
break;
}
}
const Real dlnu = (std::log(numaxd) - std::log(numind)) / n_nubinsd;
const Real nu = numind * std::exp((n + 0.5) * dlnu);
ppack_r(b, ph::energy(), n) = hd * nu;
}

if constexpr (ST == SourceType::emission) {
dejbn -= ppack_r(b, ph::weight(), n);
Expand All @@ -212,9 +281,22 @@ TaskStatus SourcePhotons(T *md, const Real t_start, const Real dt) {
typedef MeshBlockData<Real> BD;
typedef MeshData<Real> D;
typedef SourceType ST;
template TaskStatus SourcePhotons<BD, ST::thermal>(BD *md, const Real t0, const Real dt);
template TaskStatus SourcePhotons<BD, ST::emission>(BD *md, const Real t0, const Real dt);
template TaskStatus SourcePhotons<D, ST::thermal>(D *md, const Real t0, const Real dt);
template TaskStatus SourcePhotons<D, ST::emission>(D *md, const Real t0, const Real dt);
typedef FrequencyType FT;
template TaskStatus SourcePhotons<BD, ST::thermal, FT::gray>(BD *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<BD, ST::emission, FT::gray>(BD *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<D, ST::thermal, FT::gray>(D *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<D, ST::emission, FT::gray>(D *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<BD, ST::thermal, FT::multigroup>(BD *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<BD, ST::emission, FT::multigroup>(BD *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<D, ST::thermal, FT::multigroup>(D *md, const Real t0,
const Real dt);
template TaskStatus SourcePhotons<D, ST::emission, FT::multigroup>(D *md, const Real t0,
const Real dt);

} // namespace jaybenne
59 changes: 46 additions & 13 deletions src/jaybenne/transport.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ namespace jaybenne {
//----------------------------------------------------------------------------------------
//! \fn TaskStatus TransportPhotons
//! \brief
template <FrequencyType FT>
TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real dt) {
namespace fj = field::jaybenne;
namespace fjh = field::jaybenne::host;
Expand All @@ -34,11 +35,26 @@ TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real d
auto pm = md->GetParentPointer();
auto &resolved_pkgs = pm->resolved_packages;
auto &jb_pkg = pm->packages.Get("jaybenne");
auto &eos = jb_pkg->Param<EOS>("eos_d");
auto &opacity = jb_pkg->Param<Opacity>("opacity_d");
auto &scattering = jb_pkg->Param<Scattering>("scattering_d");
auto &rng_pool = jb_pkg->Param<RngPool>("rng_pool");
const Real vv = jb_pkg->Param<Real>("speed_of_light");
auto &eos = jb_pkg->template Param<EOS>("eos_d");
Opacity opacity;
Scattering scattering;
MeanOpacity mopacity;
MeanScattering mscattering;
int n_nubins = -1;
Real numin = -1.;
Real numax = -1.;
if constexpr (FT == FrequencyType::gray) {
mopacity = jb_pkg->template Param<MeanOpacity>("mopacity_d");
mscattering = jb_pkg->template Param<MeanScattering>("mscattering_d");
} else if constexpr (FT == FrequencyType::multigroup) {
opacity = jb_pkg->template Param<Opacity>("opacity_d");
scattering = jb_pkg->template Param<Scattering>("scattering_d");
n_nubins = jb_pkg->template Param<int>("n_nubins");
numin = jb_pkg->template Param<Real>("numin");
numax = jb_pkg->template Param<Real>("numax");
}
auto &rng_pool = jb_pkg->template Param<RngPool>("rng_pool");
const Real vv = jb_pkg->template Param<Real>("speed_of_light");

// Create SparsePack
static auto desc =
Expand Down Expand Up @@ -111,20 +127,28 @@ TaskStatus TransportPhotons(MeshData<Real> *md, const Real t_start, const Real d
"Particle initially outside X3 logical bnds!");

// Calculate cell bounds
const Real xl = coords.Xc<parthenon::X1DIR>(ip) - 0.5 * dx_i;
const Real xu = coords.Xc<parthenon::X1DIR>(ip) + 0.5 * dx_i;
const Real yl = coords.Xc<parthenon::X2DIR>(jp) - 0.5 * dx_j;
const Real yu = coords.Xc<parthenon::X2DIR>(jp) + 0.5 * dx_j;
const Real zl = coords.Xc<parthenon::X3DIR>(kp) - 0.5 * dx_k;
const Real zu = coords.Xc<parthenon::X3DIR>(kp) + 0.5 * dx_k;
const Real xl = coords.template Xc<parthenon::X1DIR>(ip) - 0.5 * dx_i;
const Real xu = coords.template Xc<parthenon::X1DIR>(ip) + 0.5 * dx_i;
const Real yl = coords.template Xc<parthenon::X2DIR>(jp) - 0.5 * dx_j;
const Real yu = coords.template Xc<parthenon::X2DIR>(jp) + 0.5 * dx_j;
const Real zl = coords.template Xc<parthenon::X3DIR>(kp) - 0.5 * dx_k;
const Real zu = coords.template Xc<parthenon::X3DIR>(kp) + 0.5 * dx_k;

// Extract physical quantities
const Real &rho = vmesh(b, fjh::density(), kp, jp, ip);
const Real &sie = vmesh(b, fjh::sie(), kp, jp, ip);
const Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
const Real &ff = vmesh(b, fj::fleck_factor(), kp, jp, ip);
const Real ss = scattering.TotalScatteringCoefficient(rho, temp, ee);
const Real aa = opacity.AbsorptionCoefficient(rho, temp, ee);
Real ss;
Real aa;
if constexpr (FT == FrequencyType::gray) {
// TODO: use TotalScatteringCoefficient(rho, temp), when available
ss = mscattering.RosselandMeanTotalScatteringCoefficient(rho, temp);
aa = mopacity.AbsorptionCoefficient(rho, temp);
} else if constexpr (FT == FrequencyType::multigroup) {
ss = scattering.TotalScatteringCoefficient(rho, temp, ee);
aa = opacity.AbsorptionCoefficient(rho, temp, ee);
}

// reset collision indicators
bool is_absorbed = false;
Expand Down Expand Up @@ -215,4 +239,13 @@ TaskStatus CheckCompletion(MeshData<Real> *md, const Real t_end) {
}
}

//----------------------------------------------------------------------------------------
//! template instantiations
template TaskStatus TransportPhotons<FrequencyType::gray>(MeshData<Real> *md,
const Real t_start,
const Real dt);
template TaskStatus TransportPhotons<FrequencyType::multigroup>(MeshData<Real> *md,
const Real t_start,
const Real dt);

} // namespace jaybenne
Loading
Loading