30332291031

Committed 26 Jul 2026 04:28AM UTC coverage: 77.831% (+0.02%) from 77.816%

Build # 30332291031

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1117 from wdhawkins/precursor_fixes

Precursor fixes for transient problems

Coverage Stats

112 of 134 new or added lines in 13 files covered. (83.58%)

4 existing lines in 3 files now uncovered.

26510 of 34061 relevant lines covered (77.83%)

81752699.98 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

84.59

/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h"
#include "framework/field_functions/field_function_grid_based.h"
#include "framework/materials/multi_group_xs/multi_group_xs.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/utils/hdf_utils.h"
#include "framework/object_factory.h"
#include "framework/logging/log.h"
#include "framework/runtime.h"
#include "framework/data_types/allowable_range.h"
#include "framework/utils/error.h"
#include "framework/utils/timer.h"
#include "framework/utils/caliper_scopes.h"
#include "caliper/cali.h"
#include <algorithm>
#include <iomanip>
#include <fstream>
#include <cstring>
#include <cassert>
#include <memory>
#include <stdexcept>
#include <sys/stat.h>
#include <unordered_map>
#include <functional>
#include <utility>

namespace opensn
{

InputParameters
LBSProblem::GetInputParameters()
{
  InputParameters params = Problem::GetInputParameters();

  params.ChangeExistingParamToOptional("name", "LBSProblem");

  params.AddRequiredParameter<std::shared_ptr<MeshContinuum>>("mesh", "Mesh");

  params.AddRequiredParameter<unsigned int>("num_groups",
                                            "The total number of groups within the solver");

  params.AddRequiredParameterArray("groupsets",
                                   "An array of blocks each specifying the input parameters for a "
                                   "<TT>LBSGroupset</TT>.");
  params.LinkParameterToBlock("groupsets", "LBSGroupset");

  params.AddRequiredParameterArray("xs_map",
                                   "Cross-section map from block IDs to cross-section objects.");

  params.AddOptionalParameterArray<std::shared_ptr<VolumetricSource>>(
    "volumetric_sources", {}, "An array of handles to volumetric sources.");

  params.AddOptionalParameterArray<std::shared_ptr<PointSource>>(
    "point_sources", {}, "An array of point sources.");

  params.AddOptionalParameterBlock(
    "options", ParameterBlock(), "Block of options. See <TT>OptionsBlock</TT>.");
  params.LinkParameterToBlock("options", "OptionsBlock");

  params.AddOptionalParameter("use_gpus", false, "Offload the sweep computation to GPUs.");

  return params;
}

LBSProblem::LBSProblem(const InputParameters& params)
  : Problem(params),
    num_groups_(params.GetParamValue<unsigned int>("num_groups")),
    grid_(params.GetSharedPtrParam<MeshContinuum>("mesh")),
    use_gpus_(params.GetParamValue<bool>("use_gpus"))
{
  // Check system for GPU acceleration
  if (use_gpus_)
  {
#ifdef __OPENSN_WITH_GPU__
    CheckCapableDevices();
#else
    OpenSnInvalidArgument(
      GetName() + ": GPU support was requested, but OpenSn was built without CUDA enabled.");
#endif // __OPENSN_WITH_GPU__
  }

  // Initialize options
  if (params.IsParameterValid("options"))
  {
    auto options_params = LBSProblem::GetOptionsBlock();
    options_params.AssignParameters(params.GetParam("options"));
    ParseOptions(options_params);
  }

  // Set geometry type
  geometry_type_ = grid_->GetGeometryType();
  OpenSnInvalidArgumentIf(geometry_type_ == GeometryType::INVALID,
                          GetName() + ": Invalid geometry type.");

  InitializeGroupsets(params);
  InitializeSources(params);
  InitializeXSMap(params);
  InitializeMaterials();
}

const LBSOptions&
LBSProblem::GetOptions() const
{
  return options_;
}

double
LBSProblem::GetTime() const
{
  return time_;
}

void
LBSProblem::SetTime(double time)
{
  time_ = time;
}

void
LBSProblem::SetTimeStep(double dt)
{
  OpenSnInvalidArgumentIf(dt <= 0.0, GetName() + ": dt must be greater than zero.");
  dt_ = dt;
}

double
LBSProblem::GetTimeStep() const
{
  return dt_;
}

void
LBSProblem::SetTheta(double theta)
{
  OpenSnInvalidArgumentIf(theta <= 0.0 or theta > 1.0,
                          GetName() + ": theta must be in (0.0, 1.0].");
  theta_ = theta;
}

double
LBSProblem::GetTheta() const
{
  return theta_;
}

bool
LBSProblem::IsTimeDependent() const
{
  return false;
}

void
LBSProblem::SetTimeDependentMode()
{
  OpenSnLogicalError(GetName() + ": Time-dependent mode is not supported for this problem type.");
}

void
LBSProblem::SetSteadyStateMode()
{
  // Steady-state is the default for problem types without time-dependent support.
}

GeometryType
LBSProblem::GetGeometryType() const
{
  return geometry_type_;
}

unsigned int
LBSProblem::GetNumMoments() const
{
  return num_moments_;
}

unsigned int
LBSProblem::GetMaxCellDOFCount() const
{
  return max_cell_dof_count_;
}

unsigned int
LBSProblem::GetMinCellDOFCount() const
{
  return min_cell_dof_count_;
}

bool
LBSProblem::UseGPUs() const
{
  return use_gpus_;
}

unsigned int
LBSProblem::GetNumGroups() const
{
  return num_groups_;
}

unsigned int
LBSProblem::GetScatteringOrder() const
{
  return scattering_order_;
}

unsigned int
LBSProblem::GetNumPrecursors() const
{
  return num_precursors_;
}

unsigned int
LBSProblem::GetMaxPrecursorsPerMaterial() const
{
  return max_precursors_per_material_;
}

const std::vector<LBSGroupset>&
LBSProblem::GetGroupsets() const
{
  return groupsets_;
}

LBSGroupset&
LBSProblem::GetGroupset(size_t groupset_id)
{
  return groupsets_.at(groupset_id);
}

const LBSGroupset&
LBSProblem::GetGroupset(size_t groupset_id) const
{
  return groupsets_.at(groupset_id);
}

size_t
LBSProblem::GetNumGroupsets() const
{
  return groupsets_.size();
}

void
LBSProblem::AddPointSource(std::shared_ptr<PointSource> point_source)
{
  point_sources_.push_back(point_source);
  point_sources_.back()->Initialize(*this);
}

void
LBSProblem::ClearPointSources()
{
  point_sources_.clear();
}

const std::vector<std::shared_ptr<PointSource>>&
LBSProblem::GetPointSources() const
{
  return point_sources_;
}

void
LBSProblem::AddVolumetricSource(std::shared_ptr<VolumetricSource> volumetric_source)
{
  volumetric_sources_.push_back(volumetric_source);
  volumetric_sources_.back()->Initialize(*this);
}

void
LBSProblem::ClearVolumetricSources()
{
  volumetric_sources_.clear();
}

const std::vector<std::shared_ptr<VolumetricSource>>&
LBSProblem::GetVolumetricSources() const
{
  return volumetric_sources_;
}

const BlockID2XSMap&
LBSProblem::GetBlockID2XSMap() const
{
  return block_id_to_xs_map_;
}

void
LBSProblem::SetBlockID2XSMap(const BlockID2XSMap& xs_map)
{
  const BlockID2XSMap old_xs_map = block_id_to_xs_map_;
  const size_t old_max_precursors_per_material = max_precursors_per_material_;
  const auto old_precursor_new_state = precursor_new_local_;
  const auto old_precursor_old_state = precursor_old_local_;

  block_id_to_xs_map_ = xs_map;
  InitializeMaterials();

  if (options_.use_precursors)
  {
    const size_t num_cells = grid_->local_cells.size();
    const size_t new_max_precursors_per_material = max_precursors_per_material_;
    const size_t num_precursor_dofs = num_cells * new_max_precursors_per_material;

    std::vector<double> remapped_precursors_new(num_precursor_dofs, 0.0);
    std::vector<double> remapped_precursors_old(num_precursor_dofs, 0.0);
    if (old_precursor_new_state.size() == num_cells * old_max_precursors_per_material)
    {
      for (const auto& cell : grid_->local_cells)
      {
        unsigned int old_num_precursors = 0;
        if (const auto old_xs_it = old_xs_map.find(cell.block_id); old_xs_it != old_xs_map.end())
        {
          const auto& old_xs = old_xs_it->second;
          if (old_xs->IsFissionable())
            old_num_precursors = old_xs->GetPrecursors().size();
        }

        const auto& new_xs = block_id_to_xs_map_.at(cell.block_id);
        const unsigned int new_num_precursors =
          new_xs->IsFissionable() ? new_xs->GetPrecursors().size() : 0;
        const unsigned int num_precursors_to_copy =
          std::min(old_num_precursors, new_num_precursors);

        const size_t old_base = cell.local_id * old_max_precursors_per_material;
        const size_t new_base = cell.local_id * new_max_precursors_per_material;
        for (unsigned int j = 0; j < num_precursors_to_copy; ++j)
        {
          remapped_precursors_new[new_base + j] = old_precursor_new_state[old_base + j];
          remapped_precursors_old[new_base + j] = old_precursor_old_state[old_base + j];
        }
      }
    }

    precursor_new_local_ = std::move(remapped_precursors_new);
    precursor_old_local_ = std::move(remapped_precursors_old);
  }
  else
  {
    precursor_new_local_.clear();
    precursor_old_local_.clear();
  }

  ResetGPUCarriers();
  InitializeGPUExtras();
}

std::shared_ptr<MeshContinuum>
LBSProblem::GetGrid() const
{
  return grid_;
}

const SpatialDiscretization&
LBSProblem::GetSpatialDiscretization() const
{
  return *discretization_;
}

const std::vector<UnitCellMatrices>&
LBSProblem::GetUnitCellMatrices() const
{
  return unit_cell_matrices_;
}

const std::map<uint64_t, UnitCellMatrices>&
LBSProblem::GetUnitGhostCellMatrices() const
{
  return unit_ghost_cell_matrices_;
}

const std::vector<CellLBSView>&
LBSProblem::GetCellTransportViews() const
{
  return cell_transport_views_;
}

OutflowBank&
LBSProblem::GetOutflowBank()
{
  return outflow_bank_;
}

std::vector<CellOutflowView>&
LBSProblem::GetCellOutflowViews()
{
  return cell_outflow_views_;
}

const std::vector<CellOutflowView>&
LBSProblem::GetCellOutflowViews() const
{
  return cell_outflow_views_;
}

void
LBSProblem::ConfigureOutflowStorage(const bool include_internal_faces)
{
  store_internal_outflows_ = include_internal_faces;
  cell_outflow_views_.clear();
  outflow_bank_ = OutflowBank(*grid_, num_groups_, store_internal_outflows_);
  cell_outflow_views_ = outflow_bank_.GetCellOutflowViews();
  ResetGPUCarriers();
  InitializeGPUExtras();
}

const UnknownManager&
LBSProblem::GetUnknownManager() const
{
  return flux_moments_uk_man_;
}

size_t
LBSProblem::GetLocalNodeCount() const
{
  return local_node_count_;
}

std::uint64_t
LBSProblem::GetGlobalNodeCount() const
{
  return global_node_count_;
}

std::vector<double>&
LBSProblem::GetQMomentsLocal()
{
  return q_moments_local_;
}

const std::vector<double>&
LBSProblem::GetQMomentsLocal() const
{
  return q_moments_local_;
}

const std::vector<double>&
LBSProblem::GetExtSrcMomentsLocal() const
{
  return ext_src_moments_local_;
}

void
LBSProblem::SetExtSrcMomentsFrom(const std::vector<double>& ext_src_moments)
{
  if (not phi_old_local_.empty())
    OpenSnLogicalErrorIf(ext_src_moments.size() != phi_old_local_.size(),
                         "SetExtSrcMomentsFrom size mismatch. Provided size=" +
                           std::to_string(ext_src_moments.size()) +
                           ", expected local DOFs=" + std::to_string(phi_old_local_.size()) + ".");

  if (ext_src_moments_local_.empty())
  {
    ext_src_moments_local_ = ext_src_moments;
    return;
  }

  assert(ext_src_moments.size() == ext_src_moments_local_.size() &&
         "SetExtSrcMomentsFrom size mismatch.");
  ext_src_moments_local_ = ext_src_moments;
}

std::vector<double>&
LBSProblem::GetPhiOldLocal()
{
  return phi_old_local_;
}

const std::vector<double>&
LBSProblem::GetPhiOldLocal() const
{
  return phi_old_local_;
}

std::vector<double>&
LBSProblem::GetPhiNewLocal()
{
  return phi_new_local_;
}

const std::vector<double>&
LBSProblem::GetPhiNewLocal() const
{
  return phi_new_local_;
}

std::vector<double>&
LBSProblem::GetPrecursorsNewLocal()
{
  return precursor_new_local_;
}

const std::vector<double>&
LBSProblem::GetPrecursorsNewLocal() const
{
  return precursor_new_local_;
}

std::vector<double>&
LBSProblem::GetPrecursorsOldLocal()
{
  return precursor_old_local_;
}

const std::vector<double>&
LBSProblem::GetPrecursorsOldLocal() const
{
  return precursor_old_local_;
}

SetSourceFunction
LBSProblem::GetActiveSetSourceFunction() const
{
  return active_set_source_function_;
}

void
LBSProblem::SetActiveSetSourceFunction(SetSourceFunction source_function)
{
  active_set_source_function_ = std::move(source_function);
}

std::pair<std::uint64_t, std::uint64_t>
LBSProblem::GetNumPhiIterativeUnknowns()
{
  const auto& sdm = *discretization_;
  const auto num_local_phi_dofs = sdm.GetNumLocalDOFs(flux_moments_uk_man_);
  const auto num_global_phi_dofs = sdm.GetNumGlobalDOFs(flux_moments_uk_man_);

  return {num_local_phi_dofs, num_global_phi_dofs};
}

InputParameters
LBSProblem::GetOptionsBlock()
{
  InputParameters params;

  params.SetGeneralDescription("Set options from a large list of parameters");
  params.AddOptionalParameter("max_mpi_message_size",
                              32768,
                              "The maximum MPI message size used during sweep initialization.");
  params.AddOptionalParameter(
    "restart_writes_enabled", false, "Flag that controls writing of restart dumps");
  params.AddOptionalParameter("write_delayed_psi_to_restart",
                              true,
                              "Flag that controls writing of delayed angular fluxes to restarts.");
  params.AddOptionalParameter("write_angular_flux_to_restart",
                              true,
                              "Flag that controls writing angular fluxes to restart dumps when "
                              "`save_angular_flux` is enabled.");
  params.AddOptionalParameter(
    "read_restart_path", "", "Full path for reading restart dumps including file stem.");
  params.AddOptionalParameter(
    "read_initial_condition_path",
    "",
    "Full path for reading restart data as an initial condition, including file stem.");
  params.AddOptionalParameter(
    "write_restart_path", "", "Full path for writing restart dumps including file stem.");
  params.AddOptionalParameter("write_restart_time_interval",
                              0,
                              "Time interval in seconds at which restart data is to be written.");
  params.AddOptionalParameter("use_precursors", true, "Flag for using delayed neutron precursors.");
  params.AddOptionalParameter("use_source_moments",
                              false,
                              "Flag for ignoring fixed sources and selectively using source "
                              "moments obtained elsewhere.");
  params.AddOptionalParameter(
    "save_angular_flux", false, "Flag indicating whether angular fluxes are to be stored or not.");
  params.AddOptionalParameter(
    "adjoint", false, "Flag for toggling whether the solver is in adjoint mode.");
  params.AddOptionalParameter(
    "verbose_inner_iterations",
    true,
    "Flag to control verbosity of inner iterations, including WGS and AGS iterations.");
  params.AddOptionalParameter(
    "verbose_outer_iterations", true, "Flag to control verbosity of outer iterations.");
  params.AddOptionalParameter(
    "max_ags_iterations", 100, "Maximum number of across-groupset iterations.");
  params.AddOptionalParameter("ags_tolerance", 1.0e-6, "Across-groupset iterations tolerance.");
  params.AddOptionalParameter("ags_convergence_check",
                              "l2",
                              "Type of convergence check for AGS iterations. Valid values are "
                              "`\"l2\"` and '\"pointwise\"'");
  params.AddOptionalParameter("power_default_kappa",
                              3.20435e-11,
                              "Default `kappa` value (Energy released per fission) to use for "
                              "power generation when cross sections do not have `kappa` values. "
                              "Default: 3.20435e-11 Joule (corresponding to 200 MeV per fission).");
  params.AddOptionalParameter("field_function_prefix_option",
                              "prefix",
                              "Prefix option on field function names. Default: `\"prefix\"`. Can "
                              "be `\"prefix\"` or `\"solver_name\"`. By default this option uses "
                              "the value of the `field_function_prefix` parameter. If this "
                              "parameter is not set, flux field functions will be exported as "
                              "`phi_gXXX_mYYY` where `XXX` is the zero padded 3 digit group number "
                              "and `YYY` is the zero padded 3 digit moment.");
  params.AddOptionalParameter("field_function_prefix",
                              "",
                              "Prefix to use on all field functions. Default: `\"\"`. By default "
                              "this option is empty. Ff specified, flux moments will be exported "
                              "as `prefix_phi_gXXX_mYYY` where `XXX` is the zero padded 3 digit "
                              "group number and `YYY` is the zero padded 3 digit moment. The "
                              "underscore after \"prefix\" is added automatically.");
  params.ConstrainParameterRange("ags_convergence_check",
                                 AllowableRangeList::New({"l2", "pointwise"}));
  params.ConstrainParameterRange("field_function_prefix_option",
                                 AllowableRangeList::New({"prefix", "solver_name"}));
  params.ConstrainParameterRange("max_mpi_message_size", AllowableRangeLowLimit::New(1024));
  params.ConstrainParameterRange("write_restart_time_interval", AllowableRangeLowLimit::New(0));
  params.ConstrainParameterRange("max_ags_iterations", AllowableRangeLowLimit::New(0));
  params.ConstrainParameterRange("ags_tolerance", AllowableRangeLowLimit::New(1.0e-18));
  params.ConstrainParameterRange("power_default_kappa", AllowableRangeLowLimit::New(0.0, false));

  return params;
}

InputParameters
LBSProblem::GetXSMapEntryBlock()
{
  InputParameters params;
  params.SetGeneralDescription("Set the cross-section map for the solver.");
  params.AddRequiredParameterArray("block_ids", "Mesh block IDs");
  params.AddRequiredParameter<std::shared_ptr<MultiGroupXS>>("xs", "Cross-section object");
  return params;
}

void
LBSProblem::ParseOptions(const InputParameters& input)
{
  auto params = LBSProblem::GetOptionsBlock();
  params.AssignParameters(input);
  const auto& params_at_assignment = input.GetParametersAtAssignment();
  const auto& specified_params = params_at_assignment.GetNumParameters() > 0
                                   ? params_at_assignment
                                   : static_cast<const ParameterBlock&>(input);

  using OptionSetter = std::function<void(const ParameterBlock&)>;
  const std::unordered_map<std::string, OptionSetter> option_setters = {
    {"max_mpi_message_size",
     [this](const ParameterBlock& spec) { options_.max_mpi_message_size = spec.GetValue<int>(); }},
    {"restart_writes_enabled",
     [this](const ParameterBlock& spec)
     { options_.restart.writes_enabled = spec.GetValue<bool>(); }},
    {"write_delayed_psi_to_restart",
     [this](const ParameterBlock& spec)
     { options_.restart.write_delayed_psi = spec.GetValue<bool>(); }},
    {"write_angular_flux_to_restart",
     [this](const ParameterBlock& spec)
     { options_.restart.write_angular_flux = spec.GetValue<bool>(); }},
    {"read_restart_path",
     [this](const ParameterBlock& spec)
     { options_.restart.read_path = BuildRestartPath(spec.GetValue<std::string>()); }},
    {"read_initial_condition_path",
     [this](const ParameterBlock& spec)
     {
       options_.restart.read_initial_condition_path =
         BuildRestartPath(spec.GetValue<std::string>());
     }},
    {"write_restart_path",
     [this](const ParameterBlock& spec)
     { options_.restart.write_path = BuildRestartPath(spec.GetValue<std::string>()); }},
    {"write_restart_time_interval",
     [this](const ParameterBlock& spec)
     { options_.restart.write_time_interval = std::chrono::seconds(spec.GetValue<int>()); }},
    {"use_precursors",
     [this](const ParameterBlock& spec) { options_.use_precursors = spec.GetValue<bool>(); }},
    {"use_source_moments",
     [this](const ParameterBlock& spec) { options_.use_src_moments = spec.GetValue<bool>(); }},
    {"save_angular_flux",
     [this](const ParameterBlock& spec) { options_.save_angular_flux = spec.GetValue<bool>(); }},
    {"verbose_inner_iterations",
     [this](const ParameterBlock& spec)
     { options_.verbose_inner_iterations = spec.GetValue<bool>(); }},
    {"max_ags_iterations",
     [this](const ParameterBlock& spec) { options_.max_ags_iterations = spec.GetValue<int>(); }},
    {"ags_tolerance",
     [this](const ParameterBlock& spec) { options_.ags_tolerance = spec.GetValue<double>(); }},
    {"ags_convergence_check",
     [this](const ParameterBlock& spec)
     { options_.ags_pointwise_convergence = (spec.GetValue<std::string>() == "pointwise"); }},
    {"verbose_outer_iterations",
     [this](const ParameterBlock& spec)
     { options_.verbose_outer_iterations = spec.GetValue<bool>(); }},
    {"power_default_kappa",
     [this](const ParameterBlock& spec)
     { options_.power_default_kappa = spec.GetValue<double>(); }},
    {"field_function_prefix_option",
     [this](const ParameterBlock& spec)
     { options_.field_function_prefix_option = spec.GetValue<std::string>(); }},
    {"field_function_prefix",
     [this](const ParameterBlock& spec)
     { options_.field_function_prefix = spec.GetValue<std::string>(); }},
    {"adjoint", [this](const ParameterBlock& spec) { options_.adjoint = spec.GetValue<bool>(); }},
  };

  for (const auto& spec : specified_params.GetParameters())
  {
    const auto setter_it = option_setters.find(spec.GetName());
    if (setter_it != option_setters.end())
      setter_it->second(spec);
  }

  OpenSnInvalidArgumentIf(options_.restart.write_time_interval > std::chrono::seconds(0) and
                            not options_.restart.writes_enabled,
                          GetName() + ": `write_restart_time_interval>0` requires "
                                      "`restart_writes_enabled=true`.");

  OpenSnInvalidArgumentIf(options_.restart.write_time_interval > std::chrono::seconds(0) and
                            options_.restart.write_time_interval < std::chrono::seconds(30),
                          GetName() + ": `write_restart_time_interval` must be 0 (disabled) "
                                      "or at least 30 seconds.");

  OpenSnInvalidArgumentIf(options_.restart.writes_enabled and options_.restart.write_path.empty(),
                          GetName() + ": `restart_writes_enabled=true` requires a non-empty "
                                      "`write_restart_path`.");

  OpenSnInvalidArgumentIf(not options_.field_function_prefix.empty() and
                            options_.field_function_prefix_option != "prefix",
                          GetName() + ": non-empty `field_function_prefix` requires "
                                      "`field_function_prefix_option=\"prefix\"`.");

  if (options_.restart.writes_enabled)
  {
    const auto dir = options_.restart.write_path.parent_path();

    // Create restart directory if necessary.
    // If dir is empty, write path resolves relative to the working directory.
    if ((not dir.empty()) and opensn::mpi_comm.rank() == 0)
    {
      if (not std::filesystem::exists(dir))
      {
        OpenSnLogicalErrorIf(not std::filesystem::create_directories(dir),
                             GetName() + ": Failed to create restart directory " + dir.string());
      }
      else
        OpenSnLogicalErrorIf(not std::filesystem::is_directory(dir),
                             GetName() + ": Restart path exists but is not a directory " +
                               dir.string());
    }
    opensn::mpi_comm.barrier();
    options_.restart.MarkWriteComplete();
  }
}

std::filesystem::path
LBSProblem::BuildRestartPath(const std::string& path_stem)
{
  if (path_stem.empty())
    return {};

  auto path = std::filesystem::path(path_stem);
  path += std::to_string(opensn::mpi_comm.rank()) + ".restart.h5";
  return path;
}

bool
LBSProblem::ReadProblemRestartData(hid_t /*file_id*/,
                                   bool /*allow_transient_initialization_from_steady*/)
{
  return true;
}

bool
LBSProblem::WriteProblemRestartData(hid_t /*file_id*/) const
{
  return true;
}

void
LBSProblem::BuildRuntime()
{
  PrintSimHeader();
  mpi_comm.barrier();

  InitializeRuntimeCore();
  ValidateRuntimeModeConfiguration();
  InitializeSources();
}

void
LBSProblem::InitializeRuntimeCore()
{
  InitializeSpatialDiscretization();
  InitializeParrays();
  InitializeBoundaries();
  InitializeGPUExtras();
}

void
LBSProblem::ValidateRuntimeModeConfiguration() const
{
  if (options_.adjoint)
    if (IsTimeDependent())
      OpenSnInvalidArgument(GetName() + ": Time-dependent adjoint problems are not supported.");
}

void
LBSProblem::InitializeSources()
{
  // Initialize point sources
  for (auto& point_source : point_sources_)
    point_source->Initialize(*this);

  // Initialize volumetric sources
  for (auto& volumetric_source : volumetric_sources_)
    volumetric_source->Initialize(*this);
}

void
LBSProblem::PrintSimHeader()
{
  if (opensn::mpi_comm.rank() == 0)
  {
    std::stringstream outstr;
    outstr << "\n"
           << "Initializing " << GetName() << "\n\n"
           << "Scattering order    : " << scattering_order_ << "\n"
           << "Number of moments   : " << num_moments_ << "\n"
           << "Number of groups    : " << num_groups_ << "\n"
           << "Number of groupsets : " << groupsets_.size() << "\n\n";

    for (const auto& groupset : groupsets_)
    {
      outstr << "***** Groupset " << groupset.id << " *****\n"
             << "Groups:\n";
      const auto n_gs_groups = groupset.GetNumGroups();
      constexpr int groups_per_line = 12;
      for (size_t i = 0; i < n_gs_groups; ++i)
      {
        outstr << std::setw(5) << groupset.first_group + i << ' ';
        if ((i + 1) % groups_per_line == 0)
          outstr << '\n';
      }
      if (n_gs_groups > 0 && n_gs_groups % groups_per_line != 0)
        outstr << '\n';
    }

    log.Log() << outstr.str() << '\n';
  }
}

void
LBSProblem::InitializeSources(const InputParameters& params)
{
  if (params.Has("volumetric_sources"))
  {
    const auto& vol_srcs = params.GetParam("volumetric_sources");
    vol_srcs.RequireBlockTypeIs(ParameterBlockType::ARRAY);
    for (const auto& src : vol_srcs)
      volumetric_sources_.push_back(src.GetValue<std::shared_ptr<VolumetricSource>>());
  }

  if (params.Has("point_sources"))
  {
    const auto& pt_srcs = params.GetParam("point_sources");
    pt_srcs.RequireBlockTypeIs(ParameterBlockType::ARRAY);
    for (const auto& src : pt_srcs)
      point_sources_.push_back(src.GetValue<std::shared_ptr<PointSource>>());
  }
}

void
LBSProblem::InitializeGroupsets(const InputParameters& params)
{
  // Initialize groups
  OpenSnInvalidArgumentIf(num_groups_ == 0, GetName() + ": Number of groups must be > 0.");

  // Initialize groupsets
  const auto& groupsets_array = params.GetParam("groupsets");
  const size_t num_gs = groupsets_array.GetNumParameters();
  OpenSnInvalidArgumentIf(num_gs == 0, GetName() + ": At least one groupset must be specified.");
  for (size_t gs = 0; gs < num_gs; ++gs)
  {
    const auto& groupset_params = groupsets_array.GetParam(gs);
    InputParameters gs_input_params = LBSGroupset::GetInputParameters();
    gs_input_params.SetObjectType("LBSProblem:LBSGroupset");
    gs_input_params.AssignParameters(groupset_params);
    groupsets_.emplace_back(gs_input_params, gs, *this);
    if (groupsets_.back().GetNumGroups() == 0)
    {
      std::stringstream oss;
      oss << GetName() << ": No groups added to groupset " << groupsets_.back().id;
      OpenSnInvalidArgument(oss.str());
    }

    if (groupsets_.back().last_group >= num_groups_)
    {
      std::stringstream oss;
      oss << GetName() << ": Groupset " << groupsets_.back().id << " has last group "
          << groupsets_.back().last_group << ", but the problem only has " << num_groups_
          << " groups.";
      OpenSnInvalidArgument(oss.str());
    }

    if (gs > 0)
    {
      const auto& previous_groupset = groupsets_[gs - 1];
      const auto& current_groupset = groupsets_.back();
      const auto expected_first_group = previous_groupset.last_group + 1;
      if (current_groupset.first_group != expected_first_group)
      {
        std::stringstream oss;
        oss << GetName() << ": Groupset " << current_groupset.id << " starts at group "
            << current_groupset.first_group << ", but it should start at group "
            << expected_first_group << " to be consecutive with groupset " << previous_groupset.id
            << ".";
        OpenSnInvalidArgument(oss.str());
      }
    }
  }
}

void
LBSProblem::InitializeXSMap(const InputParameters& params)
{
  // Build XS map
  const auto& xs_array = params.GetParam("xs_map");
  const size_t num_xs = xs_array.GetNumParameters();
  for (size_t i = 0; i < num_xs; ++i)
  {
    const auto& item_params = xs_array.GetParam(i);
    InputParameters xs_entry_pars = GetXSMapEntryBlock();
    xs_entry_pars.AssignParameters(item_params);

    const auto& block_ids_param = xs_entry_pars.GetParam("block_ids");
    block_ids_param.RequireBlockTypeIs(ParameterBlockType::ARRAY);
    const auto& block_ids = block_ids_param.GetVectorValue<unsigned int>();
    auto xs = xs_entry_pars.GetSharedPtrParam<MultiGroupXS>("xs");
    for (const auto& block_id : block_ids)
      block_id_to_xs_map_[block_id] = xs;
  }
}

void
LBSProblem::InitializeMaterials()
{
  CaliperPhaseScope cali_setup_phase("Setup", CaliperSetupPhaseDepth());
  CALI_CXX_MARK_SCOPE("Materials");

  log.Log0Verbose1() << "Initializing Materials";

  // Create set of material ids locally relevant
  int invalid_mat_cell_count = 0;
  std::set<unsigned int> unique_block_ids;
  for (auto& cell : grid_->local_cells)
  {
    unique_block_ids.insert(cell.block_id);
    if (cell.block_id == std::numeric_limits<unsigned int>::max() or
        (block_id_to_xs_map_.find(cell.block_id) == block_id_to_xs_map_.end()))
      ++invalid_mat_cell_count;
  }
  const auto& ghost_cell_ids = grid_->cells.GetGhostGlobalIDs();
  for (uint64_t cell_id : ghost_cell_ids)
  {
    const auto& cell = grid_->cells[cell_id];
    unique_block_ids.insert(cell.block_id);
    if (cell.block_id == std::numeric_limits<unsigned int>::max() or
        (block_id_to_xs_map_.find(cell.block_id) == block_id_to_xs_map_.end()))
      ++invalid_mat_cell_count;
  }
  OpenSnLogicalErrorIf(invalid_mat_cell_count > 0,
                       std::to_string(invalid_mat_cell_count) +
                         " cells encountered with an invalid material id.");

  // Get ready for processing
  for (const auto& [blk_id, mat] : block_id_to_xs_map_)
  {
    mat->SetAdjointMode(options_.adjoint);

    OpenSnLogicalErrorIf(mat->GetNumGroups() < num_groups_,
                         "Cross-sections for block \"" + std::to_string(blk_id) +
                           "\" have fewer groups (" + std::to_string(mat->GetNumGroups()) +
                           ") than the simulation (" + std::to_string(num_groups_) + "). " +
                           "Cross-sections must have at least as many groups as the simulation.");
  }

  // Initialize precursor properties
  num_precursors_ = 0;
  max_precursors_per_material_ = 0;
  for (const auto& mat_id_xs : block_id_to_xs_map_)
  {
    const auto& xs = mat_id_xs.second;
    if (xs->IsFissionable())
    {
      num_precursors_ += xs->GetPrecursors().size();
      max_precursors_per_material_ = std::max(static_cast<unsigned int>(xs->GetPrecursors().size()),
                                              max_precursors_per_material_);
    }
  }

  const bool has_fissionable_precursors =
    std::any_of(block_id_to_xs_map_.begin(),
                block_id_to_xs_map_.end(),
                [](const auto& mat_id_xs)
                {
                  const auto& xs = mat_id_xs.second;
                  return xs->IsFissionable() and not xs->GetPrecursors().empty();
                });
  const bool has_fissionable_material =
    std::any_of(block_id_to_xs_map_.begin(),
                block_id_to_xs_map_.end(),
                [](const auto& mat_id_xs) { return mat_id_xs.second->IsFissionable(); });

  const bool has_any_precursor_data =
    std::any_of(block_id_to_xs_map_.begin(),
                block_id_to_xs_map_.end(),
                [](const auto& mat_id_xs)
                {
                  const auto& xs = mat_id_xs.second;
                  return xs->IsFissionable() and not xs->GetPrecursors().empty();
                });

  if (options_.use_precursors and has_fissionable_material and not has_any_precursor_data)
  {
    log.Log0Warning() << GetName()
                      << ": options.use_precursors is enabled, but no precursor data was found "
                         "in the active cross-section map. Running without delayed-neutron "
                         "precursor coupling.";
  }

  // check compatibility when at least one fissionable material has delayed-neutron data
  if (options_.use_precursors and has_fissionable_precursors)
  {
    for (const auto& [mat_id, xs] : block_id_to_xs_map_)
    {
      OpenSnInvalidArgumentIf(xs->IsFissionable() and xs->GetPrecursors().empty(),
                              GetName() + ": incompatible cross-section data for material id " +
                                std::to_string(mat_id) +
                                ". When options.use_precursors=true and "
                                "delayed-neutron precursor data is present for one fissionable "
                                "material, it must be present for all fissionable materials.");
    }
  }

  // Update transport views if available
  if (grid_->local_cells.size() == cell_transport_views_.size())
    for (const auto& cell : grid_->local_cells)
    {
      const auto& xs_ptr = block_id_to_xs_map_[cell.block_id];
      auto& transport_view = cell_transport_views_[cell.local_id];
      transport_view.ReassignXS(*xs_ptr);
    }

  mpi_comm.barrier();
}

void
LBSProblem::InitializeSpatialDiscretization()
{
  CALI_CXX_MARK_SCOPE("SpatialDiscretization");

  OpenSnLogicalErrorIf(not discretization_,
                       GetName() + ": Missing spatial discretization. Construct the problem "
                                   "through its factory Create(...) entry point.");
  log.Log() << program_timer.GetTimeString() << " Initializing spatial discretization metadata.\n";

  ComputeUnitIntegrals();
}

void
LBSProblem::ComputeUnitIntegrals()
{
  CALI_CXX_MARK_SCOPE("UnitIntegrals");

  log.Log() << program_timer.GetTimeString() << " Computing unit integrals.\n";
  const auto& sdm = *discretization_;

  const size_t num_local_cells = grid_->local_cells.size();
  unit_cell_matrices_.resize(num_local_cells);

  for (const auto& cell : grid_->local_cells)
    unit_cell_matrices_[cell.local_id] =
      ComputeUnitCellIntegrals(sdm, cell, grid_->GetCoordinateSystem());

  const auto ghost_ids = grid_->cells.GetGhostGlobalIDs();
  for (auto ghost_id : ghost_ids)
    unit_ghost_cell_matrices_[ghost_id] =
      ComputeUnitCellIntegrals(sdm, grid_->cells[ghost_id], grid_->GetCoordinateSystem());

  // Assessing global unit cell matrix storage
  std::array<size_t, 2> num_local_ucms = {unit_cell_matrices_.size(),
                                          unit_ghost_cell_matrices_.size()};
  std::array<size_t, 2> num_global_ucms = {0, 0};

  mpi_comm.all_reduce(num_local_ucms.data(), 2, num_global_ucms.data(), mpi::op::sum<size_t>());

  opensn::mpi_comm.barrier();
  log.Log() << program_timer.GetTimeString() << " Ghost cell unit cell-matrix ratio: "
            << static_cast<double>(num_global_ucms[1]) * 100 /
                 static_cast<double>(num_global_ucms[0])
            << "%";
  log.Log() << program_timer.GetTimeString() << " Cell matrices computed.";
}

void
LBSProblem::InitializeParrays()
{
  CALI_CXX_MARK_SCOPE("ParallelArrays");

  log.Log() << program_timer.GetTimeString() << " Initializing parallel arrays."
            << " G=" << num_groups_ << " M=" << num_moments_ << std::endl;

  // Initialize unknown
  // structure
  flux_moments_uk_man_.unknowns.clear();
  for (unsigned int m = 0; m < num_moments_; ++m)
  {
    flux_moments_uk_man_.AddUnknown(UnknownType::VECTOR_N, num_groups_);
    flux_moments_uk_man_.unknowns.back().name = "m" + std::to_string(m);
  }

  // Compute local # of dof
  local_node_count_ = discretization_->GetNumLocalNodes();
  global_node_count_ = discretization_->GetNumGlobalNodes();

  // Compute num of unknowns
  size_t local_unknown_count = local_node_count_ * num_groups_ * num_moments_;
  log.LogAllVerbose1() << "LBS Number of phi unknowns: " << local_unknown_count;

  // Size local vectors
  q_moments_local_.assign(local_unknown_count, 0.0);
  phi_old_local_.assign(local_unknown_count, 0.0);
  phi_new_local_.assign(local_unknown_count, 0.0);

  // Setup precursor vector
  if (options_.use_precursors)
  {
    size_t num_precursor_dofs = grid_->local_cells.size() * max_precursors_per_material_;
    precursor_new_local_.assign(num_precursor_dofs, 0.0);
    precursor_old_local_.assign(num_precursor_dofs, 0.0);
  }

  // Initialize cell transport metadata and outflow tallies.
  size_t block_MG_counter = 0; // Counts the strides of moment and group
  min_cell_dof_count_ = std::numeric_limits<unsigned int>::max();
  max_cell_dof_count_ = 0;
  cell_transport_views_.clear();
  cell_transport_views_.reserve(grid_->local_cells.size());
  for (auto& cell : grid_->local_cells)
  {
    size_t num_nodes = discretization_->GetCellNumNodes(cell);

    // compute cell volumes
    double cell_volume = 0.0;
    const auto& IntV_shapeI = unit_cell_matrices_[cell.local_id].intV_shapeI;
    for (size_t i = 0; i < num_nodes; ++i)
      cell_volume += IntV_shapeI(i);

    size_t cell_phi_address = block_MG_counter;

    const size_t num_faces = cell.faces.size();
    std::vector<bool> face_local_flags(num_faces, true);
    std::vector<int> face_locality(num_faces, opensn::mpi_comm.rank());
    std::vector<const Cell*> neighbor_cell_ptrs(num_faces, nullptr);
    int f = 0;
    for (auto& face : cell.faces)
    {
      if (not face.has_neighbor)
      {
        face_local_flags[f] = false;
        face_locality[f] = -1;
      } // if bndry
      else
      {
        const int neighbor_partition = face.GetNeighborPartitionID(grid_.get());
        face_local_flags[f] = (neighbor_partition == opensn::mpi_comm.rank());
        face_locality[f] = neighbor_partition;
        neighbor_cell_ptrs[f] = &grid_->cells[face.neighbor_id];
      }

      ++f;
    }

    max_cell_dof_count_ = std::max(max_cell_dof_count_, static_cast<unsigned int>(num_nodes));
    min_cell_dof_count_ = std::min(min_cell_dof_count_, static_cast<unsigned int>(num_nodes));
    cell_transport_views_.emplace_back(cell_phi_address,
                                       num_nodes,
                                       num_groups_,
                                       num_moments_,
                                       *block_id_to_xs_map_[cell.block_id],
                                       cell_volume,
                                       face_local_flags,
                                       face_locality,
                                       neighbor_cell_ptrs);
    block_MG_counter += num_nodes * num_groups_ * num_moments_;
  } // for local cell
  cell_outflow_views_.clear();
  outflow_bank_ = OutflowBank(*grid_, num_groups_, store_internal_outflows_);
  cell_outflow_views_ = outflow_bank_.GetCellOutflowViews();

  // Populate grid nodal mappings
  // This is used in the Flux Data Structure (FLUDS).
  grid_nodal_mappings_.clear();
  grid_nodal_mappings_.reserve(grid_->local_cells.size());
  for (auto& cell : grid_->local_cells)
  {
    CellFaceNodalMapping cell_nodal_mapping;
    cell_nodal_mapping.reserve(cell.faces.size());

    for (auto& face : cell.faces)
    {
      std::vector<short> face_node_mapping;
      std::vector<short> cell_node_mapping;
      int adj_face_idx = -1;

      if (face.has_neighbor)
      {
        grid_->FindAssociatedVertices(face, face_node_mapping);
        grid_->FindAssociatedCellVertices(face, cell_node_mapping);
        adj_face_idx = face.GetNeighborAdjacentFaceIndex(grid_.get());
      }

      cell_nodal_mapping.emplace_back(adj_face_idx, face_node_mapping, cell_node_mapping);
    }

    grid_nodal_mappings_.push_back(cell_nodal_mapping);
  }

  // Get grid localized communicator set
  grid_local_comm_set_ = grid_->MakeMPILocalCommunicatorSet();

  opensn::mpi_comm.barrier();
  log.Log() << program_timer.GetTimeString() << " Done with parallel arrays." << std::endl;
}

#ifndef __OPENSN_WITH_GPU__
void
LBSProblem::InitializeGPUExtras()
{
}

void
LBSProblem::ResetGPUCarriers()
{
}

void
LBSProblem::CheckCapableDevices()
{
}
#endif // __OPENSN_WITH_GPU__

std::vector<double>
LBSProblem::MakeSourceMomentsFromPhi()
{
  CALI_CXX_MARK_SCOPE("Source/MomentsFromPhi");

  auto num_local_dofs = discretization_->GetNumLocalDOFs(flux_moments_uk_man_);

  std::vector<double> source_moments(num_local_dofs, 0.0);
  for (auto& groupset : groupsets_)
  {
    active_set_source_function_(groupset,
                                source_moments,
                                phi_new_local_,
                                APPLY_AGS_SCATTER_SOURCES | APPLY_WGS_SCATTER_SOURCES |
                                  APPLY_AGS_FISSION_SOURCES | APPLY_WGS_FISSION_SOURCES |
                                  APPLY_PREVIOUS_PRECURSOR_SOURCES);
  }

  return source_moments;
}

LBSProblem::~LBSProblem()
{
  ResetGPUCarriers();
}

void
LBSProblem::ZeroPhi()
{
  std::fill(phi_old_local_.begin(), phi_old_local_.end(), 0.0);
  std::fill(phi_new_local_.begin(), phi_new_local_.end(), 0.0);
}

void
LBSProblem::CopyPhiNewToOld()
{
  assert(phi_old_local_.size() == phi_new_local_.size() && "Phi vectors size mismatch.");
  phi_old_local_ = phi_new_local_;
}

void
LBSProblem::SetPhiOldFrom(const std::vector<double>& phi_old)
{
  assert(phi_old.size() == phi_old_local_.size() && "SetPhiOldFrom size mismatch.");
  phi_old_local_ = phi_old;
}

void
LBSProblem::SetPrecursorsOldFrom(const std::vector<double>& precursors_old)
{
  assert(precursors_old.size() == precursor_old_local_.size() &&
         "SetPrecursorsOldFrom size mismatch.");
  precursor_old_local_ = precursors_old;
}

void
LBSProblem::SetPhiNewFrom(const std::vector<double>& phi_new)
{
  assert(phi_new.size() == phi_new_local_.size() && "SetPhiNewFrom size mismatch.");
  phi_new_local_ = phi_new;
}

void
LBSProblem::ScalePhiOld(double factor)
{
  for (auto& value : phi_old_local_)
    value *= factor;
}

void
LBSProblem::ScalePhiNew(double factor)
{
  for (auto& value : phi_new_local_)
    value *= factor;
}

void
LBSProblem::ZeroQMoments()
{
  assert(q_moments_local_.size() == phi_old_local_.size() && "Q moments/Phi size mismatch.");
  std::fill(q_moments_local_.begin(), q_moments_local_.end(), 0.0);
}

void
LBSProblem::ScaleQMoments(double factor)
{
  for (auto& value : q_moments_local_)
    value *= factor;
}

void
LBSProblem::SetQMomentsFrom(const std::vector<double>& q_moments)
{
  assert(q_moments.size() == q_moments_local_.size() && "SetQMomentsFrom size mismatch.");
  q_moments_local_ = q_moments;
}

void
LBSProblem::ScalePrecursors(double factor)
{
  for (auto& value : precursor_new_local_)
    value *= factor;
}

void
LBSProblem::ZeroPrecursors()
{
  std::fill(precursor_new_local_.begin(), precursor_new_local_.end(), 0.0);
}

void
LBSProblem::ZeroExtSrcMoments()
{
  std::fill(ext_src_moments_local_.begin(), ext_src_moments_local_.end(), 0.0);
}

void
LBSProblem::ScaleExtSrcMoments(double factor)
{
  for (auto& value : ext_src_moments_local_)
    value *= factor;
}

void
LBSProblem::SetAdjoint(bool adjoint)
{
  if (adjoint)
    if (IsTimeDependent())
      OpenSnInvalidArgument(GetName() + ": Time-dependent adjoint problems are not supported.");

  const bool mode_changed = (adjoint != options_.adjoint);
  if (not mode_changed)
    return;

  options_.adjoint = adjoint;

  // Reinitialize materials to obtain the proper forward/adjoint cross sections.
  InitializeMaterials();

  // Forward and adjoint sources are fundamentally different.
  point_sources_.clear();
  volumetric_sources_.clear();
  ClearBoundaries();

  // Reset all solution vectors.
  ZeroPhi();
  ResetDerivedSolutionVectors();
  ZeroPrecursors();
}

void
LBSProblem::SetForward()
{
  SetAdjoint(false);
}

bool
LBSProblem::IsAdjoint() const
{
  return options_.adjoint;
}

} // namespace opensn

Open-Sn / opensn / 30332291031

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous