Implementation of nonlinear diffusion acceleration

benchmarks/picca_2016/multi_slab/multi_slab.material_map

@@ -0,0 +1 @@
+0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0

benchmarks/picca_2016/multi_slab/multi_slab.prm

@@ -0,0 +1,21 @@
+set problem dimension                        = 1
+set transport model                          = saaf
+set number of groups                         = 1
+set angular quadrature name                  = gauss_legendre
+set angular quadrature order                 = 4
+set do eigenvalue calculations               = true
+set reflective boundary names                =
+
+set x, y, z max values of boundary locations = 21
+set number of cells for x, y, z directions   = 21
+set uniform refinements                      = 4
+set number of materials                      = 2
+
+set finite element polynomial degree         = 1
+
+set output file name base                    = multi_slab
+
+subsection material ID map
+set material id file name map                = 0: material_a, 1: material_b
+set material id file name                    = multi_slab.material_map
+end

benchmarks/picca_2016/multi_slab/multi_slab_rayleigh.prm

@@ -0,0 +1,22 @@
+set problem dimension                        = 1
+set transport model                          = saaf
+set number of groups                         = 1
+set angular quadrature name                  = gauss_legendre
+set angular quadrature order                 = 4
+set do eigenvalue calculations               = true
+set reflective boundary names                =
+
+set x, y, z max values of boundary locations = 21
+set number of cells for x, y, z directions   = 21
+set uniform refinements                      = 4
+set number of materials                      = 2
+
+set finite element polynomial degree         = 1
+set k_effective updater type                 = rayleigh quotient
+
+set output file name base                    = multi_slab_rayleigh
+
+subsection material ID map
+set material id file name map                = 0: material_a, 1: material_b
+set material id file name                    = multi_slab.material_map
+end

benchmarks/picca_2016/multi_slab/multi_slab_rayleigh_nda.prm

@@ -0,0 +1,23 @@
+set problem dimension                        = 1
+set transport model                          = saaf
+set number of groups                         = 1
+set angular quadrature name                  = gauss_legendre
+set angular quadrature order                 = 4
+set do eigenvalue calculations               = true
+set reflective boundary names                =
+
+set x, y, z max values of boundary locations = 21
+set number of cells for x, y, z directions   = 21
+set uniform refinements                      = 4
+set number of materials                      = 2
+
+set finite element polynomial degree         = 1
+set k_effective updater type                 = rayleigh quotient
+
+set output file name base                    = multi_slab_rayleigh
+set do nda                                   = true
+
+subsection material ID map
+set material id file name map                = 0: material_a, 1: material_b
+set material id file name                    = multi_slab.material_map
+end

codecov.yml

@@ -1,3 +1,11 @@
+coverage:
+ status:
+  project:
+   default:
+    target: 95%
+  patch:
+   default:
+    target: 95%
 ignore:
  - "**/tests/*"
  - "**/BART/inc/**/*"

src/domain/definition.h

@@ -13,7 +13,7 @@
 #include "domain/definition_i.h"
 #include "domain/finite_element/finite_element_i.h"
 #include "domain/mesh/mesh_i.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "system/system_types.h"
 
 namespace bart {

src/domain/definition_i.h

@@ -9,7 +9,7 @@
 #include <deal.II/lac/vector.h>
 
 #include "domain/domain_types.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "system/system_types.h"
 #include "utility/has_description.h"
 

src/domain/finite_element/finite_element_gaussian.h

@@ -4,7 +4,7 @@
 #include <memory>
 
 #include "domain/finite_element/finite_element.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 
 namespace bart {
 

src/domain/finite_element/tests/finite_element_gaussian_test.cpp

@@ -9,7 +9,7 @@
 #include <deal.II/fe/fe_values.h>
 #include <gtest/gtest.h>
 
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "finite_element_test.hpp"
 #include "test_helpers/gmock_wrapper.h"
 

src/domain/mesh/mesh_cartesian.cpp

@@ -10,7 +10,7 @@
 #include <deal.II/base/utilities.h>
 
 #include "domain/mesh/factory.hpp"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 
 namespace bart::domain::mesh {
 

src/domain/mesh/tests/mesh_cartesian_test.cpp

@@ -9,7 +9,7 @@
 
 #include <deal.II/base/point.h>
 #include <deal.II/grid/tria.h>
-#include <problem/parameter_types.h>
+#include <problem/parameter_types.hpp>
 
 #include "test_helpers/test_helper_functions.h"
 #include "test_helpers/gmock_wrapper.h"

src/domain/tests/definition_mock.h

@@ -7,7 +7,7 @@
 #include <deal.II/lac/vector.h>
 
 #include "domain/definition_i.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 
 #include "test_helpers/gmock_wrapper.h"
 

src/domain/tests/definition_test.cpp

@@ -13,7 +13,7 @@
 #include "test_helpers/gmock_wrapper.h"
 #include "domain/mesh/tests/mesh_mock.h"
 #include "domain/finite_element/tests/finite_element_mock.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 
 namespace {
 

src/formulation/updater/diffusion_updater.hpp

@@ -12,7 +12,7 @@
 #include "formulation/updater/scattering_source_updater_i.h"
 #include "formulation/updater/fission_source_updater_i.h"
 #include "quadrature/quadrature_set_i.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "utility/has_description.h"
 
 namespace bart {

src/formulation/updater/saaf_updater.h

@@ -11,7 +11,7 @@
 #include "formulation/updater/scattering_source_updater_i.h"
 #include "formulation/updater/fission_source_updater_i.h"
 #include "quadrature/quadrature_set_i.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "system/solution/solution_types.h"
 #include "utility/has_description.h"
 

src/formulation/updater/tests/updater_tests.h

@@ -14,7 +14,7 @@
 #include "test_helpers/gmock_wrapper.h"
 #include "test_helpers/test_helper_functions.h"
 #include "test_helpers/dealii_test_domain.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 
 namespace bart {
 

src/framework/builder/framework_builder.cpp

@@ -49,6 +49,7 @@
 
 // Iteration classes
 #include "iteration/initializer/initialize_fixed_terms_once.h"
+#include "iteration/initializer/initialize_fixed_terms_reset_moments.hpp"
 #include "iteration/group/group_solve_iteration.h"
 #include "iteration/group/group_source_iteration.h"
 #include "iteration/outer/outer_power_iteration.hpp"
@@ -350,18 +351,31 @@ auto FrameworkBuilder<dim>::BuildGroupSolution(const int n_angles) -> std::uniqu
 template<int dim>
 auto FrameworkBuilder<dim>::BuildInitializer(const std::shared_ptr<FixedTermUpdater>& updater_ptr,
                                              const int total_groups,
-                                             const int total_angles) -> std::unique_ptr<Initializer> {
+                                             const int total_angles,
+                                             const InitializerName implementation) -> std::unique_ptr<Initializer> {
   ReportBuildingComponant("Initializer");
 
   std::unique_ptr<Initializer> return_ptr = nullptr;
 
-  using InitializeOnceType = iteration::initializer::InitializeFixedTermsOnce;
+  if (implementation == InitializerName::kInitializeFixedTermsOnce) {
+    using InitializeOnceType = iteration::initializer::InitializeFixedTermsOnce;
+    return_ptr = std::move(std::make_unique<InitializeOnceType>(updater_ptr, total_groups, total_angles));
+  } else if (implementation == InitializerName::kInitializeFixedTermsAndResetMoments) {
+    using InitializeAndResetMoments = iteration::initializer::InitializeFixedTermsResetMoments;
+    return_ptr = std::move(std::make_unique<InitializeAndResetMoments>(updater_ptr, total_groups, total_angles));
+  }
 
-  return_ptr = std::move(std::make_unique<InitializeOnceType>(updater_ptr, total_groups, total_angles));
   ReportBuildSuccess(return_ptr->description());
   return return_ptr;
 }
 
+template<int dim>
+auto FrameworkBuilder<dim>::BuildInitializer(const std::shared_ptr<FixedTermUpdater>& updater_ptr,
+                                             const int total_groups,
+                                             const int total_angles) -> std::unique_ptr<Initializer> {
+  return BuildInitializer(updater_ptr, total_groups, total_angles, InitializerName::kInitializeFixedTermsOnce);
+}
+
 template<int dim>
 auto FrameworkBuilder<dim>::BuildKEffectiveUpdater() -> std::unique_ptr<KEffectiveUpdater> {
   using ReturnType = eigenvalue::k_effective::UpdaterViaRayleighQuotient;

src/framework/builder/framework_builder.hpp

@@ -56,6 +56,7 @@ class FrameworkBuilder : public data_port::StatusDataPort, public FrameworkBuild
 
   using typename FrameworkBuilderI<dim>::DiffusionFormulationImpl;
   using typename FrameworkBuilderI<dim>::MomentCalculatorImpl;
+  using typename FrameworkBuilderI<dim>::InitializerName;
 
   using typename FrameworkBuilderI<dim>::ColorStatusPair;
   using typename FrameworkBuilderI<dim>::ColorStatusInstrument;
@@ -100,6 +101,10 @@ class FrameworkBuilder : public data_port::StatusDataPort, public FrameworkBuild
   [[nodiscard]] auto BuildInitializer(const std::shared_ptr<FixedTermUpdater>&,
                                       const int total_groups,
                                       const int total_angles) -> std::unique_ptr<Initializer> override;
+  [[nodiscard]] auto BuildInitializer(const std::shared_ptr<FixedTermUpdater>&,
+                                      const int total_groups,
+                                      const int total_angles,
+                                      const InitializerName) -> std::unique_ptr<Initializer> override;
   [[nodiscard]] auto BuildKEffectiveUpdater() -> std::unique_ptr<KEffectiveUpdater> override;
   [[nodiscard]] auto BuildKEffectiveUpdater(
       const std::shared_ptr<FiniteElement>&,

src/framework/builder/framework_builder_i.hpp

@@ -26,10 +26,11 @@
 #include "iteration/group/group_solve_iteration_i.h"
 #include "iteration/outer/outer_iteration_i.hpp"
 #include "iteration/subroutine/subroutine_i.hpp"
+#include "iteration/initializer/factory.hpp"
 #include "quadrature/calculators/spherical_harmonic_moments_i.h"
 #include "quadrature/calculators/angular_flux_integrator_i.hpp"
 #include "quadrature/quadrature_set_i.h"
-#include "problem/parameter_types.h"
+#include "problem/parameter_types.hpp"
 #include "solver/group/single_group_solver_i.h"
 #include "system/moments/spherical_harmonic_i.h"
 #include "system/moments/spherical_harmonic_types.h"
@@ -79,6 +80,7 @@ class FrameworkBuilderI {
 
   // Implementation specifiers
   using DiffusionFormulationImpl = formulation::DiffusionFormulationImpl;
+  using InitializerName = iteration::initializer::InitializerName;
   using MomentCalculatorImpl = quadrature::MomentCalculatorImpl;
   using SubroutineName = iteration::subroutine::SubroutineName;
 
@@ -126,6 +128,10 @@ class FrameworkBuilderI {
   virtual auto BuildInitializer(const std::shared_ptr<FixedTermUpdater>&,
                                 const int total_groups,
                                 const int total_angles) -> std::unique_ptr<Initializer> = 0;
+  virtual auto BuildInitializer(const std::shared_ptr<FixedTermUpdater>&,
+                                const int total_groups,
+                                const int total_angles,
+                                const InitializerName) -> std::unique_ptr<Initializer> = 0;
   virtual auto BuildKEffectiveUpdater() -> std::unique_ptr<KEffectiveUpdater> = 0;
   virtual auto BuildKEffectiveUpdater(
       const std::shared_ptr<FiniteElement>&,

src/framework/builder/framework_validator.cpp

@@ -32,26 +32,10 @@ auto FrameworkValidator::Parse(const framework::FrameworkParameters parameters)
   }
 
   // Check for required angular solution storage
-  if (!parameters.reflective_boundaries.empty() &&
-      parameters.equation_type == problem::EquationType::kSelfAdjointAngularFlux) {
-    needed_parts_.insert(FrameworkPart::AngularSolutionStorage);
-  }
-}
-
-void FrameworkValidator::Parse(const problem::ParametersI& to_parse) {
-  needed_parts_ = {FrameworkPart::ScatteringSourceUpdate};
-  if (to_parse.IsEigenvalueProblem())
-    needed_parts_.insert(FrameworkPart::FissionSourceUpdate);
-
-  const auto reflective_boundary = to_parse.ReflectiveBoundary();
-  bool has_reflective = std::any_of(
-      reflective_boundary.begin(),
-      reflective_boundary.end(),
-      [](std::pair<problem::Boundary, bool> pair){ return pair.second; });
-
-  if (has_reflective &&
-      to_parse.TransportModel() == problem::EquationType::kSelfAdjointAngularFlux) {
+  if (const bool equation_type_is_saaf{ parameters.equation_type == problem::EquationType::kSelfAdjointAngularFlux};
+      parameters.use_nda_ || (!parameters.reflective_boundaries.empty() && equation_type_is_saaf)) {
     needed_parts_.insert(FrameworkPart::AngularSolutionStorage);
+    if (parameters.use_nda_) AssertThrow(equation_type_is_saaf, dealii::ExcMessage("NDA requires angular solve"))
   }
 }
 

src/framework/builder/framework_validator.hpp

@@ -20,7 +20,6 @@ class FrameworkValidator : public FrameworkValidatorI,
  public:
   auto AddPart(const FrameworkPart to_add) -> FrameworkValidator& override;
   auto Parse(const framework::FrameworkParameters) -> void override;
-  auto Parse(const problem::ParametersI& to_parse) -> void override;
   auto ReportValidation() -> void override;
 
   [[nodiscard]] auto HasNeededParts() const -> bool  override { return needed_parts_.size() > 0; }

src/framework/builder/framework_validator_i.hpp

@@ -19,7 +19,6 @@ class FrameworkValidatorI {
   virtual ~FrameworkValidatorI() = default;
   virtual auto AddPart(const FrameworkPart to_add) -> FrameworkValidatorI& = 0;
   virtual auto Parse(const framework::FrameworkParameters) -> void = 0;
-  virtual auto Parse(const problem::ParametersI& to_parse) -> void = 0;
   virtual auto ReportValidation() -> void = 0;
 
   virtual auto HasNeededParts() const -> bool = 0;

src/framework/builder/tests/framework_builder_integration_tests.cpp

@@ -36,6 +36,7 @@
 #include "solver/group/single_group_solver.h"
 #include "system/solution/mpi_group_angular_solution.h"
 #include "iteration/initializer/initialize_fixed_terms_once.h"
+#include "iteration/initializer/initialize_fixed_terms_reset_moments.hpp"
 #include "iteration/group/group_source_iteration.h"
 #include "iteration/subroutine/get_scalar_flux_from_framework.hpp"
 #include "system/system_types.h"
@@ -877,9 +878,9 @@ TYPED_TEST(FrameworkBuilderIntegrationTest, BuildSystem) {
 class FrameworkBuilderIntegrationNonDimTest
  : public FrameworkBuilderIntegrationTest<bart::testing::OneD> {};
 
-TEST_F(FrameworkBuilderIntegrationNonDimTest, BuildInitializer) {
-  auto fixed_updater_ptr =
-      std::make_shared<formulation::updater::FixedUpdaterMock>();
+TEST_F(FrameworkBuilderIntegrationNonDimTest, BuildDefaultInitializer) {
+  using InitializerName = iteration::initializer::InitializerName;
+  auto fixed_updater_ptr = std::make_shared<formulation::updater::FixedUpdaterMock>();
 
   using ExpectedType = iteration::initializer::InitializeFixedTermsOnce;
   const int total_groups = bart::test_helpers::RandomDouble(1, 10);
@@ -895,4 +896,42 @@ TEST_F(FrameworkBuilderIntegrationNonDimTest, BuildInitializer) {
   EXPECT_EQ(dynamic_ptr->total_groups(), total_groups);
 }
 
+TEST_F(FrameworkBuilderIntegrationNonDimTest, BuildInitializerFixedOnceSpecified) {
+  using InitializerName = iteration::initializer::InitializerName;
+  auto fixed_updater_ptr = std::make_shared<formulation::updater::FixedUpdaterMock>();
+
+  using ExpectedType = iteration::initializer::InitializeFixedTermsOnce;
+  const int total_groups = bart::test_helpers::RandomDouble(1, 10);
+  const int total_angles = total_groups + 1;
+
+  auto initializer_ptr = this->test_builder_ptr_->BuildInitializer(fixed_updater_ptr,
+                                                                   total_groups,
+                                                                   total_angles,
+                                                                   InitializerName::kInitializeFixedTermsOnce);
+  auto dynamic_ptr = dynamic_cast<ExpectedType*>(initializer_ptr.get());
+  ASSERT_NE(initializer_ptr, nullptr);
+  ASSERT_NE(dynamic_ptr, nullptr);
+  EXPECT_EQ(dynamic_ptr->total_angles(), total_angles);
+  EXPECT_EQ(dynamic_ptr->total_groups(), total_groups);
+}
+
+TEST_F(FrameworkBuilderIntegrationNonDimTest, BuildInitializerResetMoments) {
+  using InitializerName = iteration::initializer::InitializerName;
+  auto fixed_updater_ptr = std::make_shared<formulation::updater::FixedUpdaterMock>();
+
+  using ExpectedType = iteration::initializer::InitializeFixedTermsResetMoments;
+  const int total_groups = bart::test_helpers::RandomDouble(1, 10);
+  const int total_angles = total_groups + 1;
+
+  auto initializer_ptr = this->test_builder_ptr_->BuildInitializer(fixed_updater_ptr,
+                                                                   total_groups,
+                                                                   total_angles,
+                                                                   InitializerName::kInitializeFixedTermsAndResetMoments);
+  auto dynamic_ptr = dynamic_cast<ExpectedType*>(initializer_ptr.get());
+  ASSERT_NE(initializer_ptr, nullptr);
+  ASSERT_NE(dynamic_ptr, nullptr);
+  EXPECT_EQ(dynamic_ptr->total_angles(), total_angles);
+  EXPECT_EQ(dynamic_ptr->total_groups(), total_groups);
+}
+
 } // namespace