Adding transient solver

framework/materials/multi_group_xs/multi_group_xs.cc

@@ -66,7 +66,7 @@ MultiGroupXS::Combine(
                          "All cross sections being combined must have the same group structure.");
 
     // Increment number of precursors
-    n_precs += mgxs.num_precursors_;
+    n_precs += xs->GetNumPrecursors();
   } // for cross section
 
   // Check that the fissile and precursor densities are greater than
@@ -146,7 +146,8 @@ MultiGroupXS::Combine(
       if (xsecs[x]->IsFissionable())
       {
         mgxs.sigma_f_[g] += sig_f[g];
-        mgxs.chi_[g] += ff_i * chi[g];
+        if (not chi.empty())
+          mgxs.chi_[g] += ff_i * chi[g];
         mgxs.nu_sigma_f_[g] += sig_f[g];
         for (size_t gp = 0; gp < mgxs.num_groups_; ++gp)
           mgxs.production_matrix_[g][gp] += F[g][gp];
@@ -183,11 +184,12 @@ MultiGroupXS::Combine(
     }
 
     // Set inverse velocity data
-    if (x == 0 and xsecs[x]->GetInverseVelocity().empty())
+    if (x == 0 and not xsecs[x]->GetInverseVelocity().empty())
       mgxs.inv_velocity_ = xsecs[x]->GetInverseVelocity();
-    OpenSnLogicalErrorIf(
-      xsecs[x]->GetInverseVelocity() != mgxs.inv_velocity_,
-      "All cross sections being combined must have the same group-wise velocities.");
+    if (not mgxs.inv_velocity_.empty())
+      OpenSnLogicalErrorIf(
+        xsecs[x]->GetInverseVelocity() != mgxs.inv_velocity_,
+        "All cross sections being combined must have the same group-wise velocities.");
 
     // Combine transfer matrices
 

modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc

@@ -514,6 +514,12 @@ DiscreteOrdinatesProblem::EnableTimeDependentMode()
   SetSweepChunkMode(SweepChunkMode::TimeDependent);
 }
 
+void
+DiscreteOrdinatesProblem::ReinitializeSolverSchemes()
+{
+  InitializeSolverSchemes();
+}
+
 void
 DiscreteOrdinatesProblem::InitializeBoundaries()
 {

modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h

@@ -41,6 +41,8 @@ class DiscreteOrdinatesProblem : public LBSProblem
     return SetSweepChunk(groupset);
   }
   void EnableTimeDependentMode();
+  /// Rebuild WGS/AGS solver schemes (e.g., after changing sweep chunk mode).
+  void ReinitializeSolverSchemes();
 
   /// Static registration based constructor.
   explicit DiscreteOrdinatesProblem(const InputParameters& params);

modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/sweep_wgs_context.cc

@@ -48,6 +48,18 @@ SweepWGSContext::SweepWGSContext(DiscreteOrdinatesProblem& do_problem,
 {
 }
 
+void
+SweepWGSContext::RebuildAngularFluxFromConvergedPhi(bool include_rhs_time_term)
+{
+  const auto scope = lhs_src_scope | rhs_src_scope;
+  set_source_function(groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);
+
+  sweep_chunk->IncludeRHSTimeTerm(include_rhs_time_term);
+  ApplyInverseTransportOperator(scope);
+  LBSVecOps::GSScopedCopyPrimarySTLvectors(
+    do_problem, groupset, PhiSTLOption::PHI_NEW, PhiSTLOption::PHI_OLD);
+}
+
 void
 SweepWGSContext::SetPreconditioner(KSP& solver)
 {
@@ -137,17 +149,7 @@ SweepWGSContext::PostSolveCallback()
       (groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_RICHARDSON and
        groupset.max_iterations > 1))
   {
-    const auto scope = lhs_src_scope | rhs_src_scope;
-    set_source_function(
-      groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);
-
-    // Add RHS time term (tau*psi^n)
-    if (sweep_chunk->IsTimeDependent())
-      sweep_chunk->IncludeRHSTimeTerm(true);
-
-    ApplyInverseTransportOperator(scope);
-    LBSVecOps::GSScopedCopyPrimarySTLvectors(
-      do_problem, groupset, PhiSTLOption::PHI_NEW, PhiSTLOption::PHI_OLD);
+    RebuildAngularFluxFromConvergedPhi(sweep_chunk->IsTimeDependent());
   }
 
   if (log_info)

modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/sweep_wgs_context.h

@@ -28,6 +28,7 @@ struct SweepWGSContext : public WGSContext
   void ApplyInverseTransportOperator(SourceFlags scope) override;
 
   void PostSolveCallback() override;
+  void RebuildAngularFluxFromConvergedPhi(bool include_rhs_time_term);
 
   std::shared_ptr<SweepChunk> sweep_chunk;
   SweepScheduler sweep_scheduler;

modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.cc

@@ -144,8 +144,8 @@ LBSProblem::GetTimeStep() const
 void
 LBSProblem::SetTheta(double theta)
 {
-  if (theta < 0.0 or theta > 1.0)
-    throw std::runtime_error(GetName() + " theta must be between 0.0 and 1.0.");
+  if (theta <= 0.0 or theta > 1.0)
+    throw std::runtime_error(GetName() + " theta must be in (0.0, 1.0].");
   theta_ = theta;
 }
 
@@ -414,6 +414,12 @@ LBSProblem::GetActiveSetSourceFunction() const
   return active_set_source_function_;
 }
 
+void
+LBSProblem::SetActiveSetSourceFunction(SetSourceFunction source_function)
+{
+  active_set_source_function_ = std::move(source_function);
+}
+
 std::shared_ptr<AGSLinearSolver>
 LBSProblem::GetAGSSolver()
 {

modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h

@@ -208,6 +208,7 @@ class LBSProblem : public Problem
   const std::vector<double>& GetDensitiesLocal() const;
 
   SetSourceFunction GetActiveSetSourceFunction() const;
+  void SetActiveSetSourceFunction(SetSourceFunction source_function);
 
   std::shared_ptr<AGSLinearSolver> GetAGSSolver();
 

modules/linear_boltzmann_solvers/lbs_problem/source_functions/transient_source_function.cc

@@ -2,14 +2,13 @@
 // SPDX-License-Identifier: MIT
 
 #include "modules/linear_boltzmann_solvers/lbs_problem/source_functions/transient_source_function.h"
+#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
 
 namespace opensn
 {
 
-TransientSourceFunction::TransientSourceFunction(const LBSProblem& lbs_problem,
-                                                 double& ref_dt,
-                                                 SteppingMethod& method)
-  : SourceFunction(lbs_problem), dt_(ref_dt), method_(method)
+TransientSourceFunction::TransientSourceFunction(const LBSProblem& lbs_problem)
+  : SourceFunction(lbs_problem), lbs_problem_(lbs_problem)
 {
 }
 
@@ -19,16 +18,7 @@ TransientSourceFunction::DelayedFission(const PrecursorList& precursors,
                                         const std::vector<double>& nu_delayed_sigma_f,
                                         const double* phi) const
 {
-  const auto& BackwardEuler = SteppingMethod::IMPLICIT_EULER;
-  const auto& CrankNicolson = SteppingMethod::CRANK_NICOLSON;
-
-  double theta = 0.7;
-  if (method_ == BackwardEuler)
-    theta = 1.0;
-  else if (method_ == CrankNicolson)
-    theta = 0.5;
-
-  const double eff_dt = theta * dt_;
+  const double eff_dt = lbs_problem_.GetTheta() * lbs_problem_.GetTimeStep();
 
   double value = 0.0;
   if (apply_ags_fission_src_)

modules/linear_boltzmann_solvers/lbs_problem/source_functions/transient_source_function.h

@@ -4,7 +4,6 @@
 #pragma once
 
 #include "modules/linear_boltzmann_solvers/lbs_problem/source_functions/source_function.h"
-#include "framework/math/math_time_stepping.h"
 
 namespace opensn
 {
@@ -20,16 +19,15 @@ class TransientSourceFunction : public SourceFunction
    * Constructor for the transient source function. The only difference as compared to a steady
    * source function is the treatment of delayed fission.
    */
-  TransientSourceFunction(const LBSProblem& lbs_problem, double& ref_dt, SteppingMethod& method);
+  explicit TransientSourceFunction(const LBSProblem& lbs_problem);
 
   double DelayedFission(const PrecursorList& precursors,
                         const double& rho,
                         const std::vector<double>& nu_delayed_sigma_f,
                         const double* phi) const override;
 
 private:
-  double& dt_;
-  SteppingMethod& method_;
+  const LBSProblem& lbs_problem_;
 };
 
 } // namespace opensn

modules/linear_boltzmann_solvers/solvers/time_dependent_solver.cc

@@ -68,17 +68,16 @@ TimeDependentSourceSolver::Initialize()
 
   lbs_problem_->Initialize();
   ags_solver_ = lbs_problem_->GetAGSSolver();
+  initialized_ = true;
 }
 
 void
 TimeDependentSourceSolver::Execute()
 {
+  if (not initialized_)
+    throw std::runtime_error(GetName() + ": Initialize must be called before Execute.");
   const double t0 = current_time_;
   const double tf = stop_time_;
-  const double dt_nom = lbs_problem_->GetTimeStep();
-
-  if (dt_nom <= 0.0)
-    throw std::runtime_error(GetName() + ": dt must be positive");
 
   if (tf < t0)
     throw std::runtime_error(GetName() + ": stop_time must be >= current_time");
@@ -98,14 +97,17 @@ TimeDependentSourceSolver::Execute()
       break;
     }
 
-    double step_dt = (remaining < dt_nom) ? remaining : dt_nom;
+    if (pre_advance_callback_)
+      pre_advance_callback_();
+
+    const double dt = lbs_problem_->GetTimeStep();
+    if (dt <= 0.0)
+      throw std::runtime_error(GetName() + ": dt must be positive");
+    const double step_dt = (remaining < dt) ? remaining : dt;
 
     lbs_problem_->SetTimeStep(step_dt);
     lbs_problem_->SetTime(current_time_);
 
-    if (pre_advance_callback_)
-      pre_advance_callback_();
-
     Advance();
 
     if (post_advance_callback_)
@@ -119,6 +121,8 @@ TimeDependentSourceSolver::Execute()
 void
 TimeDependentSourceSolver::Advance()
 {
+  if (not initialized_)
+    throw std::runtime_error(GetName() + ": Initialize must be called before Advance.");
   const double dt = lbs_problem_->GetTimeStep();
   const double theta = lbs_problem_->GetTheta();
 

modules/linear_boltzmann_solvers/solvers/time_dependent_solver.h

@@ -50,6 +50,7 @@ class TimeDependentSourceSolver : public Solver
   double current_time_ = 0.0;
   unsigned int step_ = 0;
   bool verbose_ = true;
+  bool initialized_ = false;
   std::function<void()> pre_advance_callback_;
   std::function<void()> post_advance_callback_;