Drift diffusion formulation

src/calculator/drift_diffusion/drift_diffusion_vector_calculator.cpp

@@ -0,0 +1,28 @@
+#include "calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp"
+
+namespace bart::calculator::drift_diffusion {
+
+namespace  {
+std::string error_preamble{ "Error in call to DriftDiffusion: " };
+} // namespace
+
+
+template<int dim>
+auto DriftDiffusionVectorCalculator<dim>::DriftDiffusion(const double scalar_flux,
+                                                         const double integrated_angular_flux,
+                                                         const Tensor& shape_gradient,
+                                                         const double sigma_t,
+                                                         const double diffusion_coefficient) const -> Tensor {
+  AssertThrow(integrated_angular_flux >= 0, dealii::ExcMessage(error_preamble + "integrated angular flux is negative"))
+  AssertThrow(sigma_t > 0, dealii::ExcMessage(error_preamble + "possible called in void, sigma_t = 0"))
+  if (scalar_flux == 0)
+    return Tensor();
+  return (integrated_angular_flux * shape_gradient / sigma_t
+      - diffusion_coefficient * shape_gradient * scalar_flux) / scalar_flux;
+}
+
+template class DriftDiffusionVectorCalculator<1>;
+template class DriftDiffusionVectorCalculator<2>;
+template class DriftDiffusionVectorCalculator<3>;
+
+} // namespace bart::calculator::drift_diffusion

src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp

@@ -0,0 +1,22 @@
+#ifndef BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_HPP_
+#define BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_HPP_
+
+#include "calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp"
+
+namespace bart::calculator::drift_diffusion {
+
+template <int dim>
+class DriftDiffusionVectorCalculator : public DriftDiffusionVectorCalculatorI<dim> {
+ public:
+  using typename DriftDiffusionVectorCalculatorI<dim>::Tensor;
+
+  [[nodiscard]] auto DriftDiffusion(const double scalar_flux,
+                                    const double integrated_angular_flux,
+                                    const Tensor &shape_gradient,
+                                    const double sigma_t,
+                                    const double diffusion_coefficient) const -> Tensor override;
+};
+
+} // namespace bart::calculator::drift_diffusion
+
+#endif //BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_HPP_

src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp

@@ -0,0 +1,40 @@
+#ifndef BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_I_HPP_
+#define BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_I_HPP_
+
+#include <deal.II/base/tensor.h>
+
+namespace bart::calculator::drift_diffusion {
+
+/*! \brief Interface for classes that calculate the drift-diffusion vector.
+ *
+ * This class has the single job of combining the terms of the drift-diffusion vector and returning the correct value.
+ * The vector is calculated for a single degree of freedom and quadrature point within the cell. Therefore, the values
+ * of the scalar flux and integrated angular flux are both scalar values.
+ *
+ * \f[
+ * D(\phi, I(\Psi), \nabla\varphi, \sigma_t, D) = \frac{1}{\phi}\left[\frac{1}{\sigma_t}I(\Psi)\nabla\varphi
+ *                                                - D\nabla\varphi\phi\right]\;,
+ * \f]
+ *
+ * where \f$\phi\f$ is the scalar flux, \f$I(\Psi)\f$ is the angular flux integrated over angle (using quadrature)
+ * \f$\sum_m w_m\hat{\Omega}_m\hat{\Omega_m}\Psi_m\f$, \f$\nabla\varphi\f$ is the gradient of the shape function,
+ * \f$\sigma_t\f$ is the total cross-section, and \f$D\f$ is the diffusion coefficent. The resulting return vector
+ * will be of length \f$d \times 1\f$ where \f$d\f$ is the dimension of the problem.
+ *
+ * @tparam dim spatial dimension, required for return tensor and shape functions.
+ *
+ */
+ template <int dim>
+class DriftDiffusionVectorCalculatorI {
+ public:
+  using Tensor = typename dealii::Tensor<1, dim>;
+  virtual auto DriftDiffusion(const double scalar_flux,
+                              const double integrated_angular_flux,
+                              const Tensor& shape_gradient,
+                              const double sigma_t,
+                              const double diffusion_coefficient) const -> Tensor = 0;
+};
+
+} // namespace bart::calculator::drift_diffusion
+
+#endif //BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_I_HPP_

src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_mock.hpp

@@ -0,0 +1,19 @@
+#ifndef BART_SRC_CALCULATOR_DRIFT_DIFFUSION_TESTS_DRIFT_DIFFUSION_VECTOR_CALCULATOR_MOCK_HPP_
+#define BART_SRC_CALCULATOR_DRIFT_DIFFUSION_TESTS_DRIFT_DIFFUSION_VECTOR_CALCULATOR_MOCK_HPP_
+
+#include "calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp"
+#include "test_helpers/gmock_wrapper.h"
+
+namespace bart::calculator::drift_diffusion {
+
+template <int dim>
+class DriftDiffusionVectorCalculatorMock : public DriftDiffusionVectorCalculatorI<dim> {
+ public:
+  using typename DriftDiffusionVectorCalculatorI<dim>::Tensor;
+  MOCK_METHOD(Tensor, DriftDiffusion, (const double scalar_flux, const double integrated_angular_flux,
+      const Tensor& shape_gradient, const double sigma_t, const double diffusion_coefficient), (const, override));
+};
+
+} // namespace bart::calculator::drift_diffusion
+
+#endif //BART_SRC_CALCULATOR_DRIFT_DIFFUSION_TESTS_DRIFT_DIFFUSION_VECTOR_CALCULATOR_MOCK_HPP_

src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp

@@ -0,0 +1,104 @@
+#include "test_helpers/gmock_wrapper.h"
+
+#include "calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp"
+#include "test_helpers/test_assertions.hpp"
+
+namespace  {
+
+using namespace bart;
+using test_helpers::AreEqual;
+
+template <typename DimensionWrapper>
+class DriftDiffusionVectorCalculatorTest : public ::testing::Test {
+ public:
+  static constexpr int dim = DimensionWrapper::value;
+  using Tensor = typename dealii::Tensor<1, dim>;
+  using TestCalculator = typename calculator::drift_diffusion::DriftDiffusionVectorCalculator<dim>;
+
+  // Test object
+  TestCalculator test_calculator_;
+
+  // Test parameters
+  const double scalar_flux_{ 5.0 };
+  const double integrated_angular_flux_{ 10.2 };
+  const double sigma_t_{ 0.25 };
+  const double diffision_coefficient_{ 2.3 };
+  Tensor shape_gradient;
+
+  auto SetUp() -> void override;
+};
+
+template <typename DimensionWrapper>
+auto DriftDiffusionVectorCalculatorTest<DimensionWrapper>::SetUp() -> void {
+  shape_gradient[0] = 1.1;
+  if (dim > 1) {
+    shape_gradient[1] = 2.2;
+    if (dim > 2) {
+      shape_gradient[2] = 3.3;
+    }
+  }
+}
+
+TYPED_TEST_SUITE(DriftDiffusionVectorCalculatorTest, bart::testing::AllDimensions);
+
+TYPED_TEST(DriftDiffusionVectorCalculatorTest, DefaultParameters) {
+  constexpr int dim{ this->dim };
+  using Tensor = dealii::Tensor<1, dim>;
+  const auto result = this->test_calculator_.DriftDiffusion(this->scalar_flux_, this->integrated_angular_flux_,
+                                                            this->shape_gradient, this->sigma_t_,
+                                                            this->diffision_coefficient_);
+  Tensor expected_result;
+  expected_result[0] = 6.446;
+  if (dim > 1) {
+    expected_result[1] = 12.892;
+    if (dim > 2) {
+      expected_result[2] = 19.338;
+    }
+  }
+  EXPECT_TRUE(AreEqual(expected_result, result));
+}
+
+TYPED_TEST(DriftDiffusionVectorCalculatorTest, ZeroSigmaT) {
+  constexpr int dim{ this->dim };
+  using Tensor = dealii::Tensor<1, dim>;
+  const double sigma_t{ 0 };
+  EXPECT_ANY_THROW({
+  [[maybe_unused]] auto result = this->test_calculator_.DriftDiffusion(this->scalar_flux_,
+                                                                       this->integrated_angular_flux_,
+                                                                       this->shape_gradient,
+                                                                       sigma_t,
+                                                                       this->diffision_coefficient_);
+                   });
+}
+
+TYPED_TEST(DriftDiffusionVectorCalculatorTest, NegativeAngularFlux) {
+  constexpr int dim{ this->dim };
+  using Tensor = dealii::Tensor<1, dim>;
+  EXPECT_ANY_THROW({
+    [[maybe_unused]] auto result = this->test_calculator_.DriftDiffusion(this->scalar_flux_,
+                                                                         -this->integrated_angular_flux_,
+                                                                         this->shape_gradient,
+                                                                         this->sigma_t_,
+                                                                         this->diffision_coefficient_);
+                   });
+}
+
+TYPED_TEST(DriftDiffusionVectorCalculatorTest, ZeroScalarFlux) {
+  constexpr int dim{ this->dim };
+  using Tensor = dealii::Tensor<1, dim>;
+  const double scalar_flux{ 0 };
+  const auto result = this->test_calculator_.DriftDiffusion(scalar_flux, this->integrated_angular_flux_,
+                                                            this->shape_gradient, this->sigma_t_,
+                                                            this->diffision_coefficient_);
+  Tensor expected_result;
+  expected_result[0] = 0;
+  if (dim > 1) {
+    expected_result[1] = 0;
+    if (dim > 2) {
+      expected_result[2] = 0;
+    }
+  }
+  EXPECT_TRUE(AreEqual(expected_result, result));
+}
+
+} // namespace

src/domain/finite_element/finite_element.cc

@@ -44,8 +44,7 @@ bool FiniteElement<dim>::SetFace(const domain::CellPtr<dim> &to_set,
   return !cell_already_set;
 }
 template<int dim>
-std::vector<double> FiniteElement<dim>::ValueAtQuadrature(
-    const system::moments::MomentVector moment) const {
+std::vector<double> FiniteElement<dim>::ValueAtQuadrature(const system::moments::MomentVector& moment) const {
 
   std::vector<double> return_vector(n_cell_quad_pts(), 0);
 

src/domain/finite_element/finite_element.h

@@ -73,8 +73,7 @@ class FiniteElement : public FiniteElementI<dim> {
     return face_values_->normal_vector(0);
   };
 
-  std::vector<double> ValueAtQuadrature(
-      const system::moments::MomentVector moment) const override;
+  std::vector<double> ValueAtQuadrature(const system::moments::MomentVector& moment) const override;
 
   std::vector<double> ValueAtFaceQuadrature(
       const dealii::Vector<double>& values_at_dofs) const override;

src/domain/finite_element/finite_element_i.h

@@ -135,8 +135,7 @@ class FiniteElementI {
    * \param moment flux moment to get the value of.
    * \return a vector holding the value of the moment at each quadrature point.
    */
-  virtual std::vector<double> ValueAtQuadrature(
-      const system::moments::MomentVector moment) const = 0;
+  virtual std::vector<double> ValueAtQuadrature(const system::moments::MomentVector& moment) const = 0;
 
   /*! \brief Get the value of an MPI Vector at the cell face quadrature points.
    *

src/domain/finite_element/tests/finite_element_mock.h

@@ -43,7 +43,7 @@ class FiniteElementMock : public FiniteElementI<dim> {
 
   MOCK_METHOD((dealii::Tensor<1, dim>), FaceNormal, (), (const, override));
 
-  MOCK_METHOD(std::vector<double>, ValueAtQuadrature, (const system::moments::MomentVector moment), (const, override));
+  MOCK_METHOD(std::vector<double>, ValueAtQuadrature, (const system::moments::MomentVector&), (const, override));
 
   MOCK_METHOD(std::vector<double>, ValueAtFaceQuadrature,
       (const dealii::Vector<double>&), (const, override));

src/formulation/factory/formulation_factories.h

@@ -14,7 +14,7 @@
 
 // Updaters
 #include "formulation/updater/saaf_updater.h"
-#include "formulation/updater/diffusion_updater.h"
+#include "formulation/updater/diffusion_updater.hpp"
 
 // Dependencies
 #include "data/cross_sections.h"

src/formulation/scalar/drift_diffusion.cpp

@@ -0,0 +1,57 @@
+#include "formulation/scalar/drift_diffusion.hpp"
+
+namespace bart::formulation::scalar {
+
+template<int dim>
+DriftDiffusion<dim>::DriftDiffusion(std::shared_ptr<FiniteElement> finite_element_ptr,
+                                    std::shared_ptr<CrossSections> cross_sections_ptr,
+                                    std::shared_ptr<DriftDiffusionCalculator> drift_diffusion_calculator_ptr)
+    : finite_element_ptr_(finite_element_ptr),
+      cross_sections_ptr_(cross_sections_ptr),
+      drift_diffusion_calculator_ptr_(drift_diffusion_calculator_ptr) {
+  std::string function_name{"formulation::scalar::DriftDiffusion constructor"};
+  AssertPointerNotNull(finite_element_ptr_.get(), "finite_element_ptr", function_name);
+  AssertPointerNotNull(cross_sections_ptr_.get(), "cross_sections_ptr", function_name);
+  AssertPointerNotNull(drift_diffusion_calculator_ptr_.get(), "drift_diffusion_calculator_ptr", function_name);
+  cell_quadrature_points_ = finite_element_ptr->n_cell_quad_pts();
+  cell_degrees_of_freedom_ = finite_element_ptr->dofs_per_cell();
+}
+template<int dim>
+auto DriftDiffusion<dim>::FillCellDriftDiffusionTerm(Matrix& to_fill,
+                                                     const CellPtr& cell_ptr,
+                                                     system::EnergyGroup group,
+                                                     const Vector& group_scalar_flux,
+                                                     const Vector& integrated_angular_flux) const -> void {
+  std::string error_prefix{"Error in DriftDiffusion<dim>::FillCellDriftDiffusionTerm: "};
+  AssertThrow(to_fill.m() == cell_quadrature_points_, dealii::ExcMessage("matrix to fill has wrong m()"))
+  AssertThrow(to_fill.n() == cell_quadrature_points_, dealii::ExcMessage("matrix to fill has wrong n()"))
+  finite_element_ptr_->SetCell(cell_ptr);
+  const auto material_id{ cell_ptr->material_id() };
+  const double sigma_t{ cross_sections_ptr_->sigma_t.at(material_id).at(group.get()) };
+  const double diffusion_coeff{ cross_sections_ptr_->diffusion_coef.at(material_id).at(group.get()) };
+
+  const auto scalar_flux_at_q{ finite_element_ptr_->ValueAtQuadrature(group_scalar_flux) };
+  const auto integrated_angular_flux_at_q{ finite_element_ptr_->ValueAtQuadrature(integrated_angular_flux) };
+
+  for (int q = 0; q < cell_quadrature_points_; ++q) {
+    const double jacobian{ finite_element_ptr_->Jacobian(q) };
+    for (int i = 0; i < cell_degrees_of_freedom_; ++i) {
+      const auto drift_diffusion{drift_diffusion_calculator_ptr_->DriftDiffusion(
+          scalar_flux_at_q.at(q),
+          integrated_angular_flux_at_q.at(q),
+          finite_element_ptr_->ShapeGradient(i, q),
+          sigma_t,
+          diffusion_coeff)};
+      const auto drift_diffusion_term{ drift_diffusion * finite_element_ptr_->ShapeValue(i, q) };
+      for (int j = 0; j < cell_degrees_of_freedom_; ++j) {
+        to_fill(i,j) += drift_diffusion_term * this->finite_element_ptr_->ShapeGradient(j, q) * jacobian;
+      }
+    }
+  }
+}
+
+template class DriftDiffusion<1>;
+template class DriftDiffusion<2>;
+template class DriftDiffusion<3>;
+
+} // namespace bart::formulation::scalar

src/formulation/scalar/drift_diffusion.hpp

@@ -0,0 +1,48 @@
+#ifndef BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_HPP_
+#define BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_HPP_
+
+#include "formulation/scalar/drift_diffusion_i.hpp"
+
+#include <memory>
+
+#include "calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp"
+#include "data/cross_sections.h"
+#include "domain/finite_element/finite_element_i.h"
+#include "utility/has_dependencies.h"
+
+namespace bart::formulation::scalar {
+
+template <int dim>
+class DriftDiffusion : public DriftDiffusionI<dim>, public utility::HasDependencies {
+ public:
+  using CrossSections = data::CrossSections;
+  using DriftDiffusionCalculator = typename calculator::drift_diffusion::DriftDiffusionVectorCalculatorI<dim>;
+  using FiniteElement = typename domain::finite_element::FiniteElementI<dim>;
+  using typename DriftDiffusionI<dim>::Matrix;
+  using typename DriftDiffusionI<dim>::CellPtr;
+  using typename DriftDiffusionI<dim>::Vector;
+
+  DriftDiffusion(std::shared_ptr<FiniteElement>,
+                 std::shared_ptr<CrossSections>,
+                 std::shared_ptr<DriftDiffusionCalculator>);
+  auto FillCellDriftDiffusionTerm(Matrix& to_fill,
+                                  const CellPtr& cell_ptr,
+                                  system::EnergyGroup group,
+                                  const Vector& group_scalar_flux,
+                                  const Vector& integrated_angular_flux) const -> void override;
+
+  auto finite_element_ptr() const -> FiniteElement* { return finite_element_ptr_.get(); }
+  auto cross_sections_ptr() const -> CrossSections* { return cross_sections_ptr_.get(); }
+  auto drift_diffusion_calculator_ptr() const -> DriftDiffusionCalculator* {
+    return drift_diffusion_calculator_ptr_.get(); }
+ private:
+  std::shared_ptr<FiniteElement> finite_element_ptr_{ nullptr };
+  std::shared_ptr<CrossSections> cross_sections_ptr_{ nullptr };
+  std::shared_ptr<DriftDiffusionCalculator> drift_diffusion_calculator_ptr_{ nullptr };
+  int cell_quadrature_points_{ 0 };
+  int cell_degrees_of_freedom_{ 0 };
+};
+
+} // namespace bart::formulation::scalar
+
+#endif //BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_HPP_