Class DiscreteOrdinatesProblem

Defined in File discrete_ordinates_problem.h

Inheritance Relationships

Base Type

public opensn::LBSProblem (Class LBSProblem)

Derived Type

public opensn::DiscreteOrdinatesCurvilinearProblem (Class DiscreteOrdinatesCurvilinearProblem)

Class Documentation

class DiscreteOrdinatesProblem : public opensn::LBSProblem 

Base class for discrete ordinates solvers.

Subclassed by opensn::DiscreteOrdinatesCurvilinearProblem

Sweep dependency data

std::map<std::shared_ptr<AngularQuadrature>, SweepOrderGroupingInfo> quadrature_unq_so_grouping_map_: Time-absorption scale 1/(theta dt) of the current DSA operators (0 when steady).

std::map<std::shared_ptr<AngularQuadrature>, std::vector<std::shared_ptr<SPDS>>> quadrature_spds_map_: Time-absorption scale 1/(theta dt) of the current DSA operators (0 when steady).

std::map<std::shared_ptr<AngularQuadrature>, std::vector<std::unique_ptr<FLUDSCommonData>>> quadrature_fluds_commondata_map_: Time-absorption scale 1/(theta dt) of the current DSA operators (0 when steady).

const std::string sweep_type_: Time-absorption scale 1/(theta dt) of the current DSA operators (0 when steady).

double dsa_time_absorption_scale_ = 0.0: Time-absorption scale 1/(theta dt) of the current DSA operators (0 when steady).

Boundary data

BoundaryBank boundary_bank_: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

std::map<uint64_t, std::shared_ptr<SweepBoundary>> sweep_boundaries_: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

std::map<uint64_t, BoundaryDefinition> boundary_definitions_: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

std::shared_ptr<BoundaryCarrier> boundary_carrier_ = nullptr: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

bool boundary_runtime_data_initialized_ = false: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

bool has_reflecting_boundaries_ = false: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

bool has_time_dependent_boundaries_ = false: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

std::set<std::uint64_t> global_boundary_ids_: Boundary IDs of the mesh on all ranks; fixed by the mesh, so computed once.

Sweep size metadata

std::size_t max_level_size_ = 0: Max level size.

std::size_t max_angleset_size_ = 0: Max angle-set size.

unsigned int max_groupset_size_ = 0: Max group-set size.

std::vector<std::vector<double>> psi_new_local_: Max level size.

std::vector<std::vector<double>> psi_old_local_: Max level size.

std::vector<double> phi_e_new_local_: Max level size.

std::optional<SweepChunkMode> sweep_chunk_mode_: Max level size.

std::string uncollided_flux_file_: Max level size.

BlockID2XSMap uncollided_xs_map_: Cross sections the uncollided flux was computed with.

std::vector<double> uncollided_flux_moments_local_: Max level size.

std::vector<double> first_collision_source_moments_local_: Max level size.

double uncollided_source_rate_ = 0.0: Max level size.

double uncollided_outflow_rate_ = 0.0: Max level size.

bool uncollided_flux_applied_ = false: Max level size.

std::shared_ptr<AGSLinearSolver> ags_solver_: Max level size.

std::vector<std::shared_ptr<WGSContext>> wgs_contexts_: Max level size.

std::vector<std::shared_ptr<LinearSolver>> wgs_solvers_: Max level size.

Construction and transport-mode controls

inline virtual bool IsTimeDependent() const override

Switch to transient mode.

If no angular flux is stored, it is reconstructed from the current steady state with the fission source divided by reconstruction_keff, the eigenvalue the state was normalized with. Use 1.0 (as the no-argument overload does) for source-driven or critical states. The value must be finite and positive; it is not retained for later switches and does not scale the physical transient source. Collective over the problem ranks.

virtual void SetTimeDependentMode() override

Switch to transient mode.

If no angular flux is stored, it is reconstructed from the current steady state with the fission source divided by reconstruction_keff, the eigenvalue the state was normalized with. Use 1.0 (as the no-argument overload does) for source-driven or critical states. The value must be finite and positive; it is not retained for later switches and does not scale the physical transient source. Collective over the problem ranks.

void SetTimeDependentMode(double reconstruction_keff)

Switch to transient mode.

If no angular flux is stored, it is reconstructed from the current steady state with the fission source divided by reconstruction_keff, the eigenvalue the state was normalized with. Use 1.0 (as the no-argument overload does) for source-driven or critical states. The value must be finite and positive; it is not retained for later switches and does not scale the physical transient source. Collective over the problem ranks.

virtual void SetSteadyStateMode() override

Switch to transient mode.

If no angular flux is stored, it is reconstructed from the current steady state with the fission source divided by reconstruction_keff, the eigenvalue the state was normalized with. Use 1.0 (as the no-argument overload does) for source-driven or critical states. The value must be finite and positive; it is not retained for later switches and does not scale the physical transient source. Collective over the problem ranks.

virtual void SetTime(double time) override

Switch to transient mode.

If no angular flux is stored, it is reconstructed from the current steady state with the fission source divided by reconstruction_keff, the eigenvalue the state was normalized with. Use 1.0 (as the no-argument overload does) for source-driven or critical states. The value must be finite and positive; it is not retained for later switches and does not scale the physical transient source. Collective over the problem ranks.

inline double GetDSATimeAbsorptionScale() const: Returns the time-absorption scale 1/(theta dt) with which the diffusion acceleration operators (WGDSA/TGDSA) were built; 0 for steady-state operators.

void UpdateDSATimeAbsorption()

Rebuilds the WGDSA/TGDSA operators if the time absorption 1/(v theta dt) of the current mode, time step, and theta differs from the one they were built with.

Called before each transient step to pick up changes of dt or theta. Mode changes rebuild the operators through RebuildRuntimeObjects().

Problem metadata and solver access

inline const std::string &GetSweepType() const: Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

inline bool UseFusedAAO() const: Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

inline unsigned int GetFusedAAOLanes() const: Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

virtual std::pair<std::uint64_t, std::uint64_t> GetNumPhiIterativeUnknowns() override: Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

std::shared_ptr<AGSLinearSolver> GetAGSSolver(): Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

std::shared_ptr<LinearSolver> GetWGSSolver(size_t groupset_id): Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

size_t GetNumWGSSolvers(): Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

void SetWGSSolveSettings(size_t groupset_id, unsigned int max_iterations, double residual_tolerance): Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

ScopedRestore PreserveWGSSolveSettings()

Saves the within-group iteration limit and residual tolerance of every groupset.

The returned guard restores them with SetWGSSolveSettings() when destroyed. Solvers hold it while acceleration methods change these settings, so that the settings outlive neither the solve nor an exception.

WGSContext &GetWGSContext(int groupset_id): Sets the iteration limit and absolute residual tolerance of a groupset’s within-group solve, on the groupset and on its WGS solver.

Balance and output helpers

BalanceTable ComputeBalanceTable(double scaling_factor = 1.0): Returns the sweep boundaries as a read-only reference.

void ComputeBalance(double scaling_factor = 1.0): Returns the sweep boundaries as a read-only reference.

const std::map<uint64_t, std::shared_ptr<SweepBoundary>> &GetSweepBoundaries() const: Returns the sweep boundaries as a read-only reference.

const std::map<uint64_t, BoundaryDefinition> &GetBoundaryDefinitions() const: Returns the sweep boundaries as a read-only reference.

std::optional<uint64_t> FindBoundaryID(const std::string &name) const

Look up the mesh boundary ID for a user-facing boundary name.

Orthogonal 2D cylindrical (RZ) meshes are addressed as rmin, rmax, zmin, and zmax, which map to the mesh boundaries xmin, xmax, ymin, and ymax. Other meshes use their boundary names directly.

Returns:: The boundary ID, or std::nullopt if the name is not valid for this problem.

virtual void ReorientAdjointSolution() override: Reorient an adjoint solution to account for backwards streaming.

void ZeroOutflowBalanceVars(LBSGroupset &groupset): Zeroes all the outflow data-structures required to compute balance.

std::vector<std::shared_ptr<FieldFunctionGridBased>> CreateAngularFluxFieldFunctionList(const std::vector<unsigned int> &groups, const std::vector<size_t> &angles)

Create angular flux field functions for the given groups and angles.

Angles are indices into the groupset quadrature associated with each group.

Construction and setup

explicit DiscreteOrdinatesProblem(const InputParameters &params): Factory-only constructor.

virtual void InitializeSpatialDiscretization() override: Creates the piecewise-linear discontinuous discretization.

virtual void BuildRuntimeData() override: Builds the sweep, acceleration, and solver data once the base runtime data is available.

virtual void RebuildRuntimeObjects() override

Also rebuilds the angular-flux storage, source function, boundaries, acceleration, and solver schemes.

Construction-fixed data (discretization, sweep orderings, FLUDS) is kept.

void InitializeSourceFunction(): Sets the source function of the current (steady-state or time-dependent) mode.

virtual void InitializeBoundaries() override: Factory-only constructor.

inline virtual bool SupportsTimeDependentMode() const: Returns false if this geometry does not support time-dependent mode.

inline virtual bool IsCurvilinear() const: Returns true for problem types that discretize the curvilinear angular-derivative terms and therefore accept cylindrical and spherical meshes.

virtual void CheckConfigurationErrors(std::vector<std::string> &errors) const override: Discrete-ordinates rules. See LBSProblem::CheckConfigurationErrors.

void CheckCSDA(std::vector<std::string> &errors) const: CSDA rules; does nothing for problems without CSDA.

void CheckNonStreamingDirections(std::vector<std::string> &errors) const: Rejects sigma_t <= 0 where an XY quadrature has a direction with no in-plane component.

Solver and sweep state

void InitializeSolverSchemes(): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

double ComputeDSATimeAbsorptionScale() const: Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

void ReinitializeSolverSchemes(): Rebuild WGS/AGS solver schemes after runtime configuration changes.

void SetSweepChunkMode(SweepChunkMode mode): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

inline void ResetSweepChunkMode(): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

void ResetMode(SweepChunkMode target_mode, double reconstruction_keff = 1.0): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

void ReconstructAngularFluxFromSteadyState(double reconstruction_keff)

Reconstructs lagged angular fluxes from the converged steady-state scalar flux via a fixed-point iteration, holding phi/q at the converged steady-state value.

Used when transitioning into time-dependent mode. If no angular flux was stored, the reconstructed distribution is rescaled per node and group to preserve the saved scalar flux.

void InitializeBoundaryCarrier(): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

void SortAngleSetsAngleIndices(): Sort the angle indices within each angle set so that one set maps exactly to another under all reflected angle mappings.

void ResetBoundaryCarrier(): Time-absorption scale 1/(theta dt) for the current mode, time step, and theta (0 when steady).

void InitializeSweepDataStructures()

This routine initializes basic sweep datastructures that are agnostic of the number of groups and essentially the groupsets.

The routine rebuilds the data structures i) quadrature_unq_so_grouping_map_, ii) quadrature_spds_map_ and iii) quadrature_fluds_templates_map_. i) is a mapping, per quadrature, to a collection of angle-index-sets where all the angles in a particular angleset share the same sweep ordering. ii) is a mapping, per quadrature, to a collection of SPDSs where each SPDS mirrors an angle-index-set in i) iii) is again a mapping, per quadrature, to a collection of Template FLUDS where each FLUDS mirrors a SPDS in ii).

void InitFluxDataStructures(LBSGroupset &groupset): Initializes fluds_ data structures.

void ResetSweepOrderings(LBSGroupset &groupset): Clears all the sweep orderings for a groupset in preparation for another.

virtual std::shared_ptr<SweepChunk> SetSweepChunk(LBSGroupset &groupset): Sets up the sweep chunk for the given discretization method.

Restart and runtime reconfiguration

virtual bool ReadProblemRestartData(hid_t file_id, bool allow_transient_initialization_from_steady, bool repartitioned) override: Parses one boundary-condition block into definitions.

virtual bool WriteProblemRestartData(hid_t file_id) const override: Parses one boundary-condition block into definitions.

virtual void ResetDerivedSolutionVectors() override: Parses one boundary-condition block into definitions.

void ParseBoundaryDefinition(const InputParameters &params, std::map<uint64_t, BoundaryDefinition> &definitions) const: Parses one boundary-condition block into definitions.

void RebuildBoundaryRuntimeData(): Parses one boundary-condition block into definitions.

virtual std::optional<std::vector<double>> ComputeDerivedFieldFunctionData(const std::string &xs_name) const override: Parses one boundary-condition block into definitions.

Factory interface

static InputParameters GetInputParameters()

static InputParameters GetBoundaryOptionsBlock()

static std::shared_ptr<DiscreteOrdinatesProblem> Create(const ParameterBlock &params)

Public Types

enum class SweepChunkMode

Values:

enumerator DEFAULT

enumerator STEADY_STATE

enumerator TIME_DEPENDENT

Public Functions

~DiscreteOrdinatesProblem() override

std::vector<std::vector<double>> &GetPsiNewLocal()

Internal angular-flux state access.

These mutable references are used by sweep chunks, transient solvers, acceleration, restart I/O, and angular-flux I/O. User-facing Python APIs expose copies or field functions rather than mutable access to this storage. Read/write access to newest updated angular flux vector.

const std::vector<std::vector<double>> &GetPsiNewLocal() const: Read access to newest updated angular flux vector.

std::vector<std::vector<double>> &GetPsiOldLocal(): Read/write access to newest updated angular flux vector.

const std::vector<std::vector<double>> &GetPsiOldLocal() const: Read access to previous angular flux vector.

std::vector<double> &GetPhiENewLocal(): Read/write access to cellwise scalar CSDA slope moments \phi_{E,g}.

const std::vector<double> &GetPhiENewLocal() const: Read access to cellwise scalar CSDA slope moments \phi_{E,g}.

void ZeroPsi()

void ZeroPhiE()

inline bool SaveAngularFluxEnabled() const

size_t GetMaxLevelSize() const

size_t GetMaxGroupsetSize() const

size_t GetMaxAngleSetSize() const

void UpdatePsiOld(): Copy psi_new to psi_old.

void SetSaveAngularFlux(bool save)

Supported runtime discrete-ordinates reconfiguration.

These methods update the dependent sweep, boundary, acceleration, and solver data owned by the problem. Lower-level sweep chunk and solver-scheme controls remain protected implementation details.

void SetBoundaryOptions(const std::vector<InputParameters> &boundary_params, bool clear_existing)

virtual void ClearBoundaries() override: Clears all the boundary conditions from the solver.

inline BoundaryCarrier *GetBoundaryCarrier()

inline const BoundaryCarrier *GetBoundaryCarrier() const

void CopyPhiAndSrcToDevice()

void CopyPhiAndOutflowBackToHost()

void TransferDeviceBoundaryData(int groupset_id, bool host_to_device, bool force = false)

Transfer data in boundary to device or vice-versa.

Parameters:

groupset_id – Groupset ID.
host_to_device – Flag indicating the direction of transfer.
force – Force update.

inline bool HasUncollidedFlux() const

void CheckNoExternalSources(const std::string &reason, std::vector<std::string> &errors) const: Appends an error, ending with reason, for each external source: volumetric and point sources, external source moments, and incoming-flux boundaries.

inline const std::vector<double> &GetFirstCollisionSourceMoments() const

inline double GetUncollidedSourceRate() const

inline double GetUncollidedOutflowRate() const

bool RemoveUncollidedFlux()

Removes a previously applied uncollided scalar flux from the iterative state.

Returns true when flux was removed.

void ComputeFluxFromUncollided(): Adds the uncollided scalar flux to the converged collided scalar flux.

Protected Types

using SweepOrderGroupingInfo = std::pair<UniqueSOGroupings, DirIDToSOMap>

Protected Functions

virtual void PrintSimHeader() override

void InitializeFCS()