diff --git a/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.cpp b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.cpp new file mode 100644 index 000000000..2f94853fc --- /dev/null +++ b/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 +auto DriftDiffusionVectorCalculator::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 diff --git a/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp new file mode 100644 index 000000000..55adc1daf --- /dev/null +++ b/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 +class DriftDiffusionVectorCalculator : public DriftDiffusionVectorCalculatorI { + public: + using typename DriftDiffusionVectorCalculatorI::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_ diff --git a/src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp new file mode 100644 index 000000000..bfe351241 --- /dev/null +++ b/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 + +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 +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_ diff --git a/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_mock.hpp b/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_mock.hpp new file mode 100644 index 000000000..915ed11c9 --- /dev/null +++ b/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 +class DriftDiffusionVectorCalculatorMock : public DriftDiffusionVectorCalculatorI { + public: + using typename DriftDiffusionVectorCalculatorI::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_ diff --git a/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp b/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp new file mode 100644 index 000000000..2984be5ad --- /dev/null +++ b/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 +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; + + // 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 +auto DriftDiffusionVectorCalculatorTest::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 diff --git a/src/domain/finite_element/finite_element.cc b/src/domain/finite_element/finite_element.cc index 16910059a..d28df911f 100644 --- a/src/domain/finite_element/finite_element.cc +++ b/src/domain/finite_element/finite_element.cc @@ -44,8 +44,7 @@ bool FiniteElement::SetFace(const domain::CellPtr &to_set, return !cell_already_set; } template -std::vector FiniteElement::ValueAtQuadrature( - const system::moments::MomentVector moment) const { +std::vector FiniteElement::ValueAtQuadrature(const system::moments::MomentVector& moment) const { std::vector return_vector(n_cell_quad_pts(), 0); diff --git a/src/domain/finite_element/finite_element.h b/src/domain/finite_element/finite_element.h index 885de1814..4682730b8 100644 --- a/src/domain/finite_element/finite_element.h +++ b/src/domain/finite_element/finite_element.h @@ -73,8 +73,7 @@ class FiniteElement : public FiniteElementI { return face_values_->normal_vector(0); }; - std::vector ValueAtQuadrature( - const system::moments::MomentVector moment) const override; + std::vector ValueAtQuadrature(const system::moments::MomentVector& moment) const override; std::vector ValueAtFaceQuadrature( const dealii::Vector& values_at_dofs) const override; diff --git a/src/domain/finite_element/finite_element_i.h b/src/domain/finite_element/finite_element_i.h index ce298a088..19d112a84 100644 --- a/src/domain/finite_element/finite_element_i.h +++ b/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 ValueAtQuadrature( - const system::moments::MomentVector moment) const = 0; + virtual std::vector ValueAtQuadrature(const system::moments::MomentVector& moment) const = 0; /*! \brief Get the value of an MPI Vector at the cell face quadrature points. * diff --git a/src/domain/finite_element/tests/finite_element_mock.h b/src/domain/finite_element/tests/finite_element_mock.h index 17355f4f4..61a2302ff 100644 --- a/src/domain/finite_element/tests/finite_element_mock.h +++ b/src/domain/finite_element/tests/finite_element_mock.h @@ -43,7 +43,7 @@ class FiniteElementMock : public FiniteElementI { MOCK_METHOD((dealii::Tensor<1, dim>), FaceNormal, (), (const, override)); - MOCK_METHOD(std::vector, ValueAtQuadrature, (const system::moments::MomentVector moment), (const, override)); + MOCK_METHOD(std::vector, ValueAtQuadrature, (const system::moments::MomentVector&), (const, override)); MOCK_METHOD(std::vector, ValueAtFaceQuadrature, (const dealii::Vector&), (const, override)); diff --git a/src/formulation/factory/formulation_factories.h b/src/formulation/factory/formulation_factories.h index 92aa0b971..854102dbd 100644 --- a/src/formulation/factory/formulation_factories.h +++ b/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" diff --git a/src/formulation/scalar/drift_diffusion.cpp b/src/formulation/scalar/drift_diffusion.cpp new file mode 100644 index 000000000..d48aaa581 --- /dev/null +++ b/src/formulation/scalar/drift_diffusion.cpp @@ -0,0 +1,57 @@ +#include "formulation/scalar/drift_diffusion.hpp" + +namespace bart::formulation::scalar { + +template +DriftDiffusion::DriftDiffusion(std::shared_ptr finite_element_ptr, + std::shared_ptr cross_sections_ptr, + std::shared_ptr 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 +auto DriftDiffusion::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::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 diff --git a/src/formulation/scalar/drift_diffusion.hpp b/src/formulation/scalar/drift_diffusion.hpp new file mode 100644 index 000000000..92a5221a6 --- /dev/null +++ b/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 + +#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 +class DriftDiffusion : public DriftDiffusionI, public utility::HasDependencies { + public: + using CrossSections = data::CrossSections; + using DriftDiffusionCalculator = typename calculator::drift_diffusion::DriftDiffusionVectorCalculatorI; + using FiniteElement = typename domain::finite_element::FiniteElementI; + using typename DriftDiffusionI::Matrix; + using typename DriftDiffusionI::CellPtr; + using typename DriftDiffusionI::Vector; + + DriftDiffusion(std::shared_ptr, + std::shared_ptr, + std::shared_ptr); + 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 finite_element_ptr_{ nullptr }; + std::shared_ptr cross_sections_ptr_{ nullptr }; + std::shared_ptr 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_ diff --git a/src/formulation/scalar/drift_diffusion_i.hpp b/src/formulation/scalar/drift_diffusion_i.hpp new file mode 100644 index 000000000..b4f6c562f --- /dev/null +++ b/src/formulation/scalar/drift_diffusion_i.hpp @@ -0,0 +1,27 @@ +#ifndef BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_HPP_ +#define BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_HPP_ + +#include +#include + +#include "domain/domain_types.h" +#include "system/system_types.h" + +namespace bart::formulation::scalar { + +template +class DriftDiffusionI { + public: + using CellPtr = typename domain::CellPtr; + using EnergyGroup = system::EnergyGroup; + using Matrix = typename dealii::FullMatrix; + using Vector = typename dealii::Vector; + virtual ~DriftDiffusionI() = default; + virtual auto FillCellDriftDiffusionTerm(Matrix& to_fill, const CellPtr&, const system::EnergyGroup, + const Vector& group_scalar_flux, + const Vector& integrated_angular_flux) const -> void = 0; +}; + +} // namespace bart::formulation::scalar + +#endif //BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_HPP_ diff --git a/src/formulation/scalar/tests/drift_diffusion_mock.hpp b/src/formulation/scalar/tests/drift_diffusion_mock.hpp new file mode 100644 index 000000000..bff0ba481 --- /dev/null +++ b/src/formulation/scalar/tests/drift_diffusion_mock.hpp @@ -0,0 +1,23 @@ +#ifndef BART_SRC_FORMULATION_SCALAR_TESTS_DRIFT_DIFFUSION_MOCK_HPP_ +#define BART_SRC_FORMULATION_SCALAR_TESTS_DRIFT_DIFFUSION_MOCK_HPP_ + +#include "formulation/scalar/drift_diffusion_i.hpp" +#include "test_helpers/gmock_wrapper.h" + +namespace bart::formulation::scalar { + +template +class DriftDiffusionMock : public DriftDiffusionI { + public: + using typename DriftDiffusionI::CellPtr; + using typename DriftDiffusionI::EnergyGroup; + using typename DriftDiffusionI::Matrix; + using typename DriftDiffusionI::Vector; + + MOCK_METHOD(void, FillCellDriftDiffusionTerm, (Matrix& to_fill, const CellPtr&, const system::EnergyGroup, + const Vector& group_scalar_flux, const Vector& integrated_angular_flux), (const, override)); +}; + +} // namespace bart::formulation::scalar + +#endif //BART_SRC_FORMULATION_SCALAR_TESTS_DRIFT_DIFFUSION_MOCK_HPP_ diff --git a/src/formulation/scalar/tests/drift_diffusion_test.cpp b/src/formulation/scalar/tests/drift_diffusion_test.cpp new file mode 100644 index 000000000..2dacbdef3 --- /dev/null +++ b/src/formulation/scalar/tests/drift_diffusion_test.cpp @@ -0,0 +1,242 @@ +#include "formulation/scalar/drift_diffusion.hpp" + +#include +#include +#include +#include + +#include "calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_mock.hpp" +#include "data/cross_sections.h" +#include "domain/finite_element/tests/finite_element_mock.h" +#include "material/tests/mock_material.h" +#include "test_helpers/gmock_wrapper.h" +#include "test_helpers/test_assertions.hpp" +#include "test_helpers/test_helper_functions.h" + +namespace { + +using namespace bart; +using test_helpers::AreEqual; + +using ::testing::NiceMock, ::testing::Return, ::testing::DoDefault, ::testing::_; +using ::testing::Ref, ::testing::AtLeast; + +template +class DriftDiffusionFormulationTest : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; + using CellPtr = typename domain::CellPtr; + using CrossSections = data::CrossSections; + using DriftDiffusionCalculator = NiceMock>; + using DriftDiffusionFormulation = formulation::scalar::DriftDiffusion; + using FiniteElement = NiceMock>; + using Material = NiceMock; + + DriftDiffusionFormulationTest() : dof_handler_(triangulation_), finite_element_(1) {}; + + // test object + std::unique_ptr test_formulation_{ nullptr }; + + std::shared_ptr drift_diffusion_calculator_mock_ptr_{ nullptr }; + std::shared_ptr finite_element_mock_ptr_{ nullptr }; + std::shared_ptr cross_sections_ptr_{ nullptr }; + + // minimum number of dealii objects to get a valid cell pointer + CellPtr cell_ptr_; + dealii::Triangulation triangulation_; + dealii::DoFHandler dof_handler_; + dealii::FE_Q finite_element_; + + // Test paraameters + const int dofs_per_cell_{ 2 }; + const int cell_quadrature_points_{ 2 }; + const int material_id_{ 1 }; + const int energy_group_{ 0 }; + dealii::FullMatrix expected_result_; + std::unordered_map> sigma_t_{{material_id_, {1.0, 2.0}}}; + std::unordered_map> diffusion_coef_{{material_id_, {0.5, 1.0}}}; + + auto SetUp() -> void override; + [[nodiscard]] auto GetShapeGradient(const int i, const int q) const -> dealii::Tensor<1, dim>; +}; + + +template +auto DriftDiffusionFormulationTest::GetShapeGradient(const int i, const int q) const +-> dealii::Tensor<1, dim> { + dealii::Tensor<1, dim> gradient_tensor; + const double val{ (i + q + 0.5) }; + switch (dim) { + case 3: + gradient_tensor[2] = 3 * val; [[fallthrough]]; + case 2: + gradient_tensor[1] = 2 * val; [[fallthrough]]; + case 1: + gradient_tensor[0] = val; + break; + } + return gradient_tensor; +} + +template +auto DriftDiffusionFormulationTest::SetUp() -> void { + drift_diffusion_calculator_mock_ptr_ = std::make_shared(); + finite_element_mock_ptr_ = std::make_shared(); + Material mock_material; + + ON_CALL(*finite_element_mock_ptr_, dofs_per_cell()).WillByDefault(Return(dofs_per_cell_)); + ON_CALL(*finite_element_mock_ptr_, n_cell_quad_pts()).WillByDefault(Return(cell_quadrature_points_)); + + for (int q = 0; q < cell_quadrature_points_; ++q) { + ON_CALL(*finite_element_mock_ptr_, Jacobian(q)).WillByDefault(Return(3 * (q + 1))); + for (int i = 0; i < dofs_per_cell_; ++i) { + auto shape_gradient{ GetShapeGradient(i, q) }; + ON_CALL(*finite_element_mock_ptr_, ShapeValue(i,q)).WillByDefault(Return(1 + i + q)); + ON_CALL(*finite_element_mock_ptr_, ShapeGradient(i,q)).WillByDefault(Return(shape_gradient)); + } + } + + dealii::Tensor<1, dim> drift_diffusion; + switch (dim) { + case 3: + drift_diffusion[2] = 1000; [[fallthrough]]; + case 2: + drift_diffusion[1] = 100; [[fallthrough]]; + case 1: + drift_diffusion[0] = 10; + break; + } + + ON_CALL(*drift_diffusion_calculator_mock_ptr_, DriftDiffusion(_, _, _, _, _)).WillByDefault(Return(drift_diffusion)); + + ON_CALL(mock_material, GetSigT()).WillByDefault(Return(sigma_t_)); + ON_CALL(mock_material, GetDiffusionCoef()).WillByDefault(Return(diffusion_coef_)); + cross_sections_ptr_ = std::make_shared(mock_material); + + dealii::GridGenerator::hyper_cube(triangulation_, 0, 1); + dof_handler_.distribute_dofs(finite_element_); + for (auto cell = dof_handler_.begin_active(); cell != dof_handler_.end(); ++cell) { + if (cell->is_locally_owned()) { + cell_ptr_ = cell; + cell_ptr_->set_material_id(material_id_); + break; + } + } + +// std::array expected_result_q_0_values{0.5, 1.5, 2.5, 3.5, +// 1.0, 3.0, 5.0, 7.0, +// 1.5, 4.5, 7.5, 10.5, +// 2.0, 6.0, 10.0, 14.0}; +// std::array expected_result_q_1_values{3.0, 5.0, 7.0, 9.0, +// 4.5, 7.5, 10.5, 13.5, +// 6.0, 10.0, 14.0, 18.0, +// 7.5, 12.5, 17.5, 22.5}; + std::array expected_result_q_0_values{0.5, 1.5, + 1.0, 3.0}; + std::array expected_result_q_1_values{3.0, 5.0, + 4.5, 7.5}; + auto expected_result_q_0_ = dealii::FullMatrix(2, 2, expected_result_q_0_values.begin()); + expected_result_q_0_ *= 3; + auto expected_result_q_1_ = dealii::FullMatrix(2, 2, expected_result_q_1_values.begin()); + expected_result_q_1_ *= 6; + expected_result_ = expected_result_q_0_; + expected_result_.add(1, expected_result_q_1_); + if (dim == 1) { + expected_result_ *= 10; + } else if (dim == 2) { + expected_result_ *= 210; + } else { + expected_result_ *= 3210; + } + + test_formulation_ = std::move(std::make_unique(finite_element_mock_ptr_, + cross_sections_ptr_, + drift_diffusion_calculator_mock_ptr_)); +} + +TYPED_TEST_SUITE(DriftDiffusionFormulationTest, bart::testing::AllDimensions); + +TYPED_TEST(DriftDiffusionFormulationTest, ConstructorAndDependencyGetters) { + EXPECT_CALL(*this->finite_element_mock_ptr_, dofs_per_cell()).WillOnce(DoDefault()); + EXPECT_CALL(*this->finite_element_mock_ptr_, n_cell_quad_pts()).WillOnce(DoDefault()); + formulation::scalar::DriftDiffusiondim> drift_diffusion( + this->finite_element_mock_ptr_, + this->cross_sections_ptr_, + this->drift_diffusion_calculator_mock_ptr_); + EXPECT_NE(drift_diffusion.finite_element_ptr(), nullptr); + EXPECT_NE(drift_diffusion.cross_sections_ptr(), nullptr); + EXPECT_NE(drift_diffusion.drift_diffusion_calculator_ptr(), nullptr); +} + +TYPED_TEST(DriftDiffusionFormulationTest, ConstructorBadDependencies) { + const int n_dependencies{ 3 }; + for (int i = 0; i < n_dependencies; ++i) { + EXPECT_ANY_THROW({ + [[maybe_unused]] formulation::scalar::DriftDiffusiondim> drift_diffusion( + i == 0 ? this->finite_element_mock_ptr_ : nullptr, + i == 1 ? this->cross_sections_ptr_ : nullptr, + i == 2 ? this->drift_diffusion_calculator_mock_ptr_ : nullptr); + }); + } +} + +TYPED_TEST(DriftDiffusionFormulationTest, FillDriftDiffusion) { + auto& finite_element_mock = *this->finite_element_mock_ptr_; + auto& drift_diffusion_calculator_mock = *this->drift_diffusion_calculator_mock_ptr_; + const dealii::Vector integrated_angular_flux_at_dofs(this->dof_handler_.n_dofs()); + const dealii::Vector scalar_flux_at_dofs(this->dof_handler_.n_dofs()); + std::vector integrated_angular_flux_at_q{ test_helpers::RandomVector(this->cell_quadrature_points_, 1, 100) }; + std::vector scalar_flux_at_q{ test_helpers::RandomVector(this->cell_quadrature_points_, 1, 100) }; + const double sigma_t{ this->sigma_t_.at(this->material_id_).at(this->energy_group_) }; + const double diffusion_coeff{ this->diffusion_coef_.at(this->material_id_).at(this->energy_group_) }; + EXPECT_CALL(finite_element_mock, ValueAtQuadrature(Ref(integrated_angular_flux_at_dofs))) + .WillOnce(Return(integrated_angular_flux_at_q)); + EXPECT_CALL(finite_element_mock, ValueAtQuadrature(Ref(scalar_flux_at_dofs))) + .WillOnce(Return(scalar_flux_at_q)); + EXPECT_CALL(finite_element_mock, SetCell(this->cell_ptr_)); + for (int q = 0; q < this->cell_quadrature_points_; ++q) { + EXPECT_CALL(finite_element_mock, Jacobian(q)).Times(AtLeast(1)).WillRepeatedly(DoDefault()); + for (int i = 0; i < this->dofs_per_cell_; ++i) { + EXPECT_CALL(finite_element_mock, ShapeGradient(i, q)).Times(AtLeast(1)).WillRepeatedly(DoDefault()); + EXPECT_CALL(finite_element_mock, ShapeValue(i, q)).Times(AtLeast(1)).WillRepeatedly(DoDefault()); + EXPECT_CALL(drift_diffusion_calculator_mock, DriftDiffusion( + scalar_flux_at_q.at(q), integrated_angular_flux_at_q.at(q), this->GetShapeGradient(i, q), + sigma_t, diffusion_coeff)) + .Times(AtLeast(1)) + .WillRepeatedly(DoDefault()); + } + } + dealii::FullMatrix cell_matrix(this->dofs_per_cell_, this->dofs_per_cell_); + cell_matrix = 0; + this->test_formulation_->FillCellDriftDiffusionTerm(cell_matrix, + this->cell_ptr_, + system::EnergyGroup(this->energy_group_), + scalar_flux_at_dofs, + integrated_angular_flux_at_dofs); + EXPECT_TRUE(AreEqual(this->expected_result_, cell_matrix)); +} + +TYPED_TEST(DriftDiffusionFormulationTest, FillCellDriftDiffusionBadMatrixSize) { + std::vector bad_sizes{this->dofs_per_cell_ - 1, this->dofs_per_cell_ + 1}; + const dealii::Vector integrated_angular_flux_at_dofs(this->dof_handler_.n_dofs()); + const dealii::Vector scalar_flux_at_dofs(this->dof_handler_.n_dofs()); + + for (const auto bad_size : bad_sizes) { + dealii::FullMatrix matrix_bad_rows(bad_size, this->dofs_per_cell_); + dealii::FullMatrix matrix_bad_cols(this->dofs_per_cell_, bad_size); + for (auto& matrix : std::array, 2>{matrix_bad_rows, matrix_bad_cols}) { + EXPECT_ANY_THROW({ + this->test_formulation_->FillCellDriftDiffusionTerm(matrix, + this->cell_ptr_, + system::EnergyGroup(this->energy_group_), + scalar_flux_at_dofs, + integrated_angular_flux_at_dofs); + }); + } + } + +} + + + +} // namespace diff --git a/src/formulation/updater/diffusion_updater.cc b/src/formulation/updater/diffusion_updater.cpp similarity index 71% rename from src/formulation/updater/diffusion_updater.cc rename to src/formulation/updater/diffusion_updater.cpp index fc332895b..b2bef1ec9 100644 --- a/src/formulation/updater/diffusion_updater.cc +++ b/src/formulation/updater/diffusion_updater.cpp @@ -1,4 +1,4 @@ -#include "formulation/updater/diffusion_updater.h" +#include "formulation/updater/diffusion_updater.hpp" namespace bart { @@ -9,10 +9,11 @@ namespace updater { template DiffusionUpdater::DiffusionUpdater( std::unique_ptr formulation_ptr, - std::unique_ptr stamper_ptr, + std::shared_ptr stamper_ptr, std::unordered_set reflective_boundaries) - : formulation_ptr_(std::move(formulation_ptr)), - stamper_ptr_(std::move(stamper_ptr)), + : FixedUpdater(stamper_ptr), + formulation_ptr_(std::move(formulation_ptr)), + stamper_ptr_(stamper_ptr), reflective_boundaries_(reflective_boundaries) { AssertThrow(formulation_ptr_ != nullptr, dealii::ExcMessage("Error in constructor of DiffusionUpdater, " @@ -30,48 +31,34 @@ DiffusionUpdater::DiffusionUpdater( } template -void DiffusionUpdater::UpdateFixedTerms( - system::System& to_update, - system::EnergyGroup energy_group, - quadrature::QuadraturePointIndex /*index*/) { - using CellPtr = domain::CellPtr; - int group = energy_group.get(); - auto fixed_matrix_ptr = - to_update.left_hand_side_ptr_->GetFixedTermPtr({group, 0}); - auto fixed_vector_ptr = - to_update.right_hand_side_ptr_->GetFixedTermPtr({group, 0}); - auto streaming_term_function = [&](formulation::FullMatrix& cell_matrix, - const CellPtr& cell_ptr) -> void { - formulation_ptr_->FillCellStreamingTerm(cell_matrix, cell_ptr, group); - }; - auto collision_term_function = [&](formulation::FullMatrix& cell_matrix, - const CellPtr& cell_ptr) -> void { - formulation_ptr_->FillCellCollisionTerm(cell_matrix, cell_ptr, group); - }; - auto fixed_term_function = [&](formulation::Vector& cell_vector, - const CellPtr& cell_ptr) -> void { - formulation_ptr_->FillCellFixedSource(cell_vector, cell_ptr, group); - }; - auto boundary_function = [&](formulation::FullMatrix& cell_matrix, - const domain::FaceIndex face_index, +auto DiffusionUpdater::SetUpFixedFunctions(system::System& /*to_update*/, + system::EnergyGroup group, + quadrature::QuadraturePointIndex /*index*/) -> void { + const auto streaming_term_function = [&, group](formulation::FullMatrix& cell_matrix, const CellPtr& cell_ptr) -> void { + formulation_ptr_->FillCellStreamingTerm(cell_matrix, cell_ptr, group.get()); }; + const auto collision_term_function = [&, group](formulation::FullMatrix& cell_matrix, const CellPtr& cell_ptr) -> void { + formulation_ptr_->FillCellCollisionTerm(cell_matrix, cell_ptr, group.get()); }; + this->fixed_matrix_functions_.push_back(streaming_term_function); + this->fixed_matrix_functions_.push_back(collision_term_function); + + auto fixed_term_function = [&, group](formulation::Vector& cell_vector, const CellPtr& cell_ptr) -> void { + formulation_ptr_->FillCellFixedSource(cell_vector, cell_ptr, group.get()); }; + this->fixed_vector_functions_.push_back(fixed_term_function); + + auto boundary_function = [&](formulation::FullMatrix& cell_matrix, const domain::FaceIndex face_index, const CellPtr& cell_ptr) -> void { using DiffusionBoundaryType = typename formulation::scalar::DiffusionI::BoundaryType; - problem::Boundary boundary = static_cast( - cell_ptr->face(face_index.get())->boundary_id()); + problem::Boundary boundary = static_cast(cell_ptr->face(face_index.get())->boundary_id()); DiffusionBoundaryType boundary_type = DiffusionBoundaryType::kVacuum; if (reflective_boundaries_.count(boundary) == 1) boundary_type = DiffusionBoundaryType::kReflective; - formulation_ptr_->FillBoundaryTerm(cell_matrix, cell_ptr, - face_index.get(), boundary_type); + + formulation_ptr_->FillBoundaryTerm(cell_matrix, cell_ptr,face_index.get(), boundary_type); }; - *fixed_matrix_ptr = 0; - *fixed_vector_ptr = 0; - stamper_ptr_->StampMatrix(*fixed_matrix_ptr, streaming_term_function); - stamper_ptr_->StampMatrix(*fixed_matrix_ptr, collision_term_function); - stamper_ptr_->StampBoundaryMatrix(*fixed_matrix_ptr, boundary_function); - stamper_ptr_->StampVector(*fixed_vector_ptr, fixed_term_function); + this->fixed_matrix_boundary_functions_.push_back(boundary_function); } + template void DiffusionUpdater::UpdateScatteringSource( system::System &to_update, diff --git a/src/formulation/updater/diffusion_updater.h b/src/formulation/updater/diffusion_updater.hpp similarity index 77% rename from src/formulation/updater/diffusion_updater.h rename to src/formulation/updater/diffusion_updater.hpp index a59aa0bef..415eb2952 100644 --- a/src/formulation/updater/diffusion_updater.h +++ b/src/formulation/updater/diffusion_updater.hpp @@ -6,6 +6,7 @@ #include "formulation/scalar/diffusion_i.h" #include "formulation/stamper_i.h" +#include "formulation/updater/fixed_updater.hpp" #include "formulation/updater/fixed_updater_i.h" #include "formulation/updater/fixed_source_updater_i.h" #include "formulation/updater/scattering_source_updater_i.h" @@ -22,22 +23,23 @@ namespace updater { template class DiffusionUpdater - : public FixedUpdaterI, public ScatteringSourceUpdaterI, + : public FixedUpdater, public ScatteringSourceUpdaterI, public FissionSourceUpdaterI, public FixedSourceUpdaterI, public utility::HasDescription { public: + using typename FixedUpdater::CellPtr; + using typename FixedUpdater::MatrixFunction; + using typename FixedUpdater::VectorFunction; + using typename FixedUpdater::MatrixBoundaryFunction; + using typename FixedUpdater::VectorBoundaryFunction; + using DiffusionFormulationType = formulation::scalar::DiffusionI; using StamperType = formulation::StamperI; DiffusionUpdater(std::unique_ptr, - std::unique_ptr, + std::shared_ptr, std::unordered_set reflective_boundaries = {}); virtual ~DiffusionUpdater() = default; - void UpdateFixedTerms( - system::System&, - system::EnergyGroup, - quadrature::QuadraturePointIndex) override; - void UpdateScatteringSource( system::System &, system::EnergyGroup, @@ -58,9 +60,11 @@ class DiffusionUpdater DiffusionFormulationType* formulation_ptr() const { return formulation_ptr_.get(); } StamperType* stamper_ptr() const { return stamper_ptr_.get(); } + protected: + auto SetUpFixedFunctions(system::System&, system::EnergyGroup, quadrature::QuadraturePointIndex) -> void override; private: std::unique_ptr formulation_ptr_; - std::unique_ptr stamper_ptr_; + std::shared_ptr stamper_ptr_; std::unordered_set reflective_boundaries_; }; diff --git a/src/formulation/updater/drift_diffusion_updater.cpp b/src/formulation/updater/drift_diffusion_updater.cpp new file mode 100644 index 000000000..dbbcb63dd --- /dev/null +++ b/src/formulation/updater/drift_diffusion_updater.cpp @@ -0,0 +1,51 @@ +#include "drift_diffusion_updater.hpp" + +namespace bart::formulation::updater { + +template +DriftDiffusionUpdater::DriftDiffusionUpdater( + std::unique_ptr diffusion_formulation_ptr, + std::unique_ptr drift_diffusion_formulation_ptr, + std::shared_ptr stamper_ptr, + std::unique_ptr integrated_flux_calculator_ptr, + std::shared_ptr high_order_moments, + AngularFluxStorageMap& angular_flux_storage_map, + std::unordered_set reflective_boundaries) + : DiffusionUpdater(std::move(diffusion_formulation_ptr), stamper_ptr, reflective_boundaries), + angular_flux_storage_map_(angular_flux_storage_map), + high_order_moments_(high_order_moments), + drift_diffusion_formulation_ptr_(std::move(drift_diffusion_formulation_ptr)), + integrated_flux_calculator_ptr_(std::move(integrated_flux_calculator_ptr)) { + std::string call_location{ "DriftDiffusionUpdaterConstructor"}; + AssertPointerNotNull(drift_diffusion_formulation_ptr_.get(), "drift diffusion formulation", call_location); + AssertPointerNotNull(integrated_flux_calculator_ptr_.get(), "integrated flux calculator", call_location); + AssertPointerNotNull(high_order_moments_.get(), "higher order moments", call_location); +} +template +auto DriftDiffusionUpdater::SetUpFixedFunctions(system::System& system, + system::EnergyGroup energy_group, + quadrature::QuadraturePointIndex quadrature_point_index) -> void { + DiffusionUpdater::SetUpFixedFunctions(system, energy_group, quadrature_point_index); + using CellPtr = domain::CellPtr; + std::map>> group_angular_flux; + for (auto& [angular_flux_index, vector_ptr] : angular_flux_storage_map_) { + auto& [solution_energy_group, solution_angle_index] = angular_flux_index; + if (solution_energy_group == energy_group) + group_angular_flux.insert({quadrature::QuadraturePointIndex(solution_angle_index.get()), vector_ptr}); + } + auto integrated_angular_flux = integrated_flux_calculator_ptr_->Integrate(group_angular_flux); + auto scalar_flux = this->high_order_moments_->GetMoment({energy_group.get(), 0, 0}); + + const auto drift_diffusion_term_function = [=, this](formulation::FullMatrix& cell_matrix, + const CellPtr& cell_ptr) -> void { + drift_diffusion_formulation_ptr_->FillCellDriftDiffusionTerm( + cell_matrix, cell_ptr, energy_group, scalar_flux, integrated_angular_flux); + }; + this->fixed_matrix_functions_.push_back(drift_diffusion_term_function); +} + +template class DriftDiffusionUpdater<1>; +template class DriftDiffusionUpdater<2>; +template class DriftDiffusionUpdater<3>; + +} // namespace bart::formulation::updater diff --git a/src/formulation/updater/drift_diffusion_updater.hpp b/src/formulation/updater/drift_diffusion_updater.hpp new file mode 100644 index 000000000..846fb5b8f --- /dev/null +++ b/src/formulation/updater/drift_diffusion_updater.hpp @@ -0,0 +1,52 @@ +#ifndef BART_SRC_FORMULATION_UPDATER_DRIFT_DIFFUSION_UPDATER_HPP_ +#define BART_SRC_FORMULATION_UPDATER_DRIFT_DIFFUSION_UPDATER_HPP_ + +#include "diffusion_updater.hpp" +#include "formulation/scalar/drift_diffusion_i.hpp" +#include "formulation/stamper_i.h" +#include "quadrature/calculators/drift_diffusion_integrated_flux_i.hpp" +#include "system/solution/solution_types.h" +#include "system/moments/spherical_harmonic_i.h" +#include "utility/has_dependencies.h" + +namespace bart::formulation::updater { + +template +class DriftDiffusionUpdater : public DiffusionUpdater, public utility::HasDependencies { + public: + using AngularFluxStorageMap = system::solution::EnergyGroupToAngularSolutionPtrMap; + using DiffusionFormulation = typename DiffusionUpdater::DiffusionFormulationType; + using DriftDiffusionFormulation = formulation::scalar::DriftDiffusionI; + using IntegratedFluxCalculator = quadrature::calculators::DriftDiffusionIntegratedFluxI; + using HighOrderMoments = system::moments::SphericalHarmonicI; + using Stamper = typename DiffusionUpdater::StamperType; + + DriftDiffusionUpdater(std::unique_ptr, + std::unique_ptr, + std::shared_ptr, + std::unique_ptr, + std::shared_ptr, + AngularFluxStorageMap&, + std::unordered_set reflective_boundaries = {}); + virtual ~DriftDiffusionUpdater() = default; + + + auto angular_flux_storage_map() const -> AngularFluxStorageMap { + return angular_flux_storage_map_; } + auto drift_diffusion_formulation_ptr() const -> DriftDiffusionFormulation* { + return drift_diffusion_formulation_ptr_.get(); } + auto high_order_moments() const -> HighOrderMoments* { + return high_order_moments_.get(); } + auto integrated_flux_calculator_ptr() const -> IntegratedFluxCalculator* { + return integrated_flux_calculator_ptr_.get(); } + protected: + auto SetUpFixedFunctions(system::System&, system::EnergyGroup, quadrature::QuadraturePointIndex) -> void override; + AngularFluxStorageMap angular_flux_storage_map_{}; + std::shared_ptr high_order_moments_; + std::unique_ptr drift_diffusion_formulation_ptr_{ nullptr }; + std::unique_ptr integrated_flux_calculator_ptr_{ nullptr }; +}; + +} // namespace bart::formulation::updater + +#endif //BART_SRC_FORMULATION_UPDATER_DRIFT_DIFFUSION_UPDATER_HPP_ diff --git a/src/formulation/updater/fixed_updater.hpp b/src/formulation/updater/fixed_updater.hpp new file mode 100644 index 000000000..faec14751 --- /dev/null +++ b/src/formulation/updater/fixed_updater.hpp @@ -0,0 +1,59 @@ +#ifndef BART_SRC_FORMULATION_UPDATER_FIXED_UPDATER_HPP_ +#define BART_SRC_FORMULATION_UPDATER_FIXED_UPDATER_HPP_ + +#include "formulation/stamper_i.h" +#include "formulation/updater/fixed_updater_i.h" + +namespace bart::formulation::updater { + +template +class FixedUpdater : public FixedUpdaterI { + public: + using Matrix = dealii::FullMatrix; + using Vector = dealii::Vector; + using CellPtr = domain::CellPtr; + using FaceIndex = domain::FaceIndex; + using MatrixFunction = std::function; + using VectorFunction = std::function; + using MatrixBoundaryFunction = std::function; + using VectorBoundaryFunction = std::function; + using Stamper = formulation::StamperI; + + FixedUpdater(std::shared_ptr stamper_ptr) : stamper_ptr_(stamper_ptr) {} + virtual ~FixedUpdater() = default; + auto UpdateFixedTerms(system::System& to_update, system::EnergyGroup group, quadrature::QuadraturePointIndex index) + -> void override { + SetUpFixedFunctions(to_update, group, index); + const int energy_group{ group.get() }; + auto fixed_matrix_ptr = to_update.left_hand_side_ptr_->GetFixedTermPtr({energy_group, 0}); + auto fixed_vector_ptr = to_update.right_hand_side_ptr_->GetFixedTermPtr({energy_group, 0}); + *fixed_matrix_ptr = 0; + *fixed_vector_ptr = 0; + + for (auto& matrix_function : fixed_matrix_functions_) + stamper_ptr_->StampMatrix(*fixed_matrix_ptr, matrix_function); + for (auto& vector_function : fixed_vector_functions_) + stamper_ptr_->StampVector(*fixed_vector_ptr, vector_function); + for (auto& matrix_boundary_function : fixed_matrix_boundary_functions_) + stamper_ptr_->StampBoundaryMatrix(*fixed_matrix_ptr, matrix_boundary_function); + // LCOV_EXCL_START + // Excluded from coverage because no current formulations utilize this. + for (auto& vector_boundary_function : fixed_vector_boundary_functions_) + stamper_ptr_->StampBoundaryVector(*fixed_vector_ptr, vector_boundary_function); + // LCOV_EXCL_STOP + } + + protected: + virtual auto SetUpFixedFunctions(system::System&, system::EnergyGroup, quadrature::QuadraturePointIndex) -> void = 0; + std::vector fixed_matrix_functions_{}; + std::vector fixed_vector_functions_{}; + std::vector fixed_matrix_boundary_functions_{}; + std::vector fixed_vector_boundary_functions_{}; + + std::shared_ptr stamper_ptr_{ nullptr }; + +}; + +} // namespace bart::formulation::updater + +#endif //BART_SRC_FORMULATION_UPDATER_FIXED_UPDATER_HPP_ diff --git a/src/formulation/updater/tests/diffusion_updater_test.cc b/src/formulation/updater/tests/diffusion_updater_test.cpp similarity index 99% rename from src/formulation/updater/tests/diffusion_updater_test.cc rename to src/formulation/updater/tests/diffusion_updater_test.cpp index 826b8ae13..b502c8c96 100644 --- a/src/formulation/updater/tests/diffusion_updater_test.cc +++ b/src/formulation/updater/tests/diffusion_updater_test.cpp @@ -1,4 +1,4 @@ -#include "formulation/updater/diffusion_updater.h" +#include "formulation/updater/diffusion_updater.hpp" #include "formulation/scalar/tests/diffusion_mock.h" #include "formulation/tests/stamper_mock.h" diff --git a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp new file mode 100644 index 000000000..74b7c5a0b --- /dev/null +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -0,0 +1,202 @@ +#include "formulation/updater/drift_diffusion_updater.hpp" + +#include + +#include "quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp" +#include "formulation/scalar/tests/diffusion_mock.h" +#include "formulation/scalar/tests/drift_diffusion_mock.hpp" +#include "formulation/tests/stamper_mock.h" +#include "formulation/updater/tests/updater_tests.h" +#include "test_helpers/gmock_wrapper.h" +#include "test_helpers/test_helper_functions.h" +#include "test_helpers/test_assertions.hpp" +#include "system/solution/solution_types.h" +#include "system/moments/tests/spherical_harmonic_mock.h" + +namespace { + +using namespace bart; +template +using UpdaterTest = bart::formulation::updater::test_helpers::UpdaterTests; + + +using ::testing::ContainerEq, ::testing::DoDefault, ::testing::_, ::testing::Ref, ::testing::Return, ::testing::ReturnRef; + +template +class FormulationUpdaterDriftDiffusionTest : public UpdaterTest { + public: + static constexpr int dim{ DimensionWrapper::value }; + using AngularFluxStorageMap = system::solution::EnergyGroupToAngularSolutionPtrMap; + using DiffusionFormulation = formulation::scalar::DiffusionMock; + using DriftDiffusionFormulation = formulation::scalar::DriftDiffusionMock; + using IntegratedFluxCalculator = quadrature::calculators::DriftDiffusionIntegratedFluxMock; + using Stamper = formulation::StamperMock; + using Boundary = problem::Boundary; + using Vector = dealii::Vector; + using GroupAngularFluxMap = std::map>; + using HighOrderMoments = system::moments::SphericalHarmonicMock; + + // Test object + using TestUpdater = formulation::updater::DriftDiffusionUpdater; + std::unique_ptr test_updater_ptr_{ nullptr }; + + // Supporting objects + AngularFluxStorageMap angular_flux_storage_map_{}; // All angular fluxes + GroupAngularFluxMap group_angular_flux_{}; // Angular fluxes for a single group (the test group) + dealii::Vector integrated_angular_flux_; + dealii::Vector group_scalar_flux_; + + // Mock observation pointers + DiffusionFormulation* diffusion_formulation_obs_ptr_{ nullptr }; + DriftDiffusionFormulation* drift_diffusion_formulation_obs_ptr_{ nullptr }; + IntegratedFluxCalculator* integrated_flux_calculator_obs_ptr_{ nullptr }; + Stamper* stamper_obs_ptr_{ nullptr }; + std::shared_ptr high_order_moments_ptr_{ nullptr }; + + // Test parameters (each should be a different value to ensure the correct values are being passed + const Vector::size_type angular_flux_size{ static_cast(test_helpers::RandomInt(5, 10)) }; + + std::unordered_set reflective_boundaries{Boundary::kXMin, Boundary::kYMax}; + + void SetUp() override; +}; + +template +void FormulationUpdaterDriftDiffusionTest::SetUp() { + UpdaterTest::SetUp(); + auto diffusion_formulation_ptr = std::make_unique(); + diffusion_formulation_obs_ptr_ = diffusion_formulation_ptr.get(); + auto drift_diffusion_formulation_ptr = std::make_unique(); + drift_diffusion_formulation_obs_ptr_ = drift_diffusion_formulation_ptr.get(); + auto integrated_flux_calculator_ptr = std::make_unique(); + integrated_flux_calculator_obs_ptr_ = integrated_flux_calculator_ptr.get(); + auto stamper_ptr = std::shared_ptr(this->MakeStamper()); + stamper_obs_ptr_ = stamper_ptr.get(); + high_order_moments_ptr_ = std::make_shared(); + + integrated_angular_flux_ = dealii::Vector(angular_flux_size); + group_scalar_flux_ = dealii::Vector(angular_flux_size); + + + using Index = system::SolutionIndex; + for (int group = 0; group < this->total_groups; ++group) { + for (int angle = 0; angle < this->total_angles; ++angle) { + auto angular_flux = std::make_shared(angular_flux_size); + Index solution_index{ system::EnergyGroup(group), system::AngleIdx(angle) }; + angular_flux_storage_map_.insert({solution_index, angular_flux}); + if (group == this->group_number) + group_angular_flux_.insert({quadrature::QuadraturePointIndex(angle), angular_flux}); + } + } + test_updater_ptr_ = std::make_unique(std::move(diffusion_formulation_ptr), + std::move(drift_diffusion_formulation_ptr), + stamper_ptr, + std::move(integrated_flux_calculator_ptr), + high_order_moments_ptr_, + angular_flux_storage_map_, + reflective_boundaries); +} + +TYPED_TEST_SUITE(FormulationUpdaterDriftDiffusionTest, bart::testing::AllDimensions); + +// Constructor should not throw, dependencies should have been saved and getters return non-null ptrs +TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorDependencyGetters) { + constexpr int dim{ this->dim }; + using DiffusionFormulation = formulation::scalar::DiffusionMock; + using DriftDiffusionFormulation = formulation::scalar::DriftDiffusionMock; + using IntegratedFluxCalculator = quadrature::calculators::DriftDiffusionIntegratedFluxMock; + using Updater = formulation::updater::DriftDiffusionUpdater; + using Stamper = formulation::StamperMock; + + std::unique_ptr test_updater{ nullptr }; + EXPECT_NO_THROW({ + test_updater = std::make_unique(std::make_unique(), + std::make_unique(), + std::make_shared(), + std::make_unique(), + this->high_order_moments_ptr_, + this->angular_flux_storage_map_); + }); + ASSERT_NE(test_updater->formulation_ptr(), nullptr); + ASSERT_NE(test_updater->drift_diffusion_formulation_ptr(), nullptr); + ASSERT_NE(test_updater->stamper_ptr(), nullptr); + ASSERT_NE(test_updater->integrated_flux_calculator_ptr(), nullptr); + EXPECT_EQ(test_updater->high_order_moments(), this->high_order_moments_ptr_.get()); + ASSERT_THAT(test_updater->angular_flux_storage_map(), ContainerEq(this->angular_flux_storage_map_)); +} + +TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorBadDependencies) { + constexpr int dim{ this->dim }; + using DiffusionFormulation = formulation::scalar::DiffusionMock; + using DriftDiffusionFormulation = formulation::scalar::DriftDiffusionMock; + using IntegratedFluxCalculator = quadrature::calculators::DriftDiffusionIntegratedFluxMock; + using Updater = formulation::updater::DriftDiffusionUpdater; + using Stamper = formulation::StamperMock; + using HighOrderMoments = system::moments::SphericalHarmonicMock; + + + std::unique_ptr test_updater{ nullptr }; + const int n_dependencies_to_check{ 5 }; + for (int i = 0; i < n_dependencies_to_check; ++i) { + EXPECT_ANY_THROW({ + test_updater = std::make_unique(i == 0 ? nullptr : std::make_unique(), + i == 1 ? nullptr : std::make_unique(), + i == 2 ? nullptr : std::make_unique(), + i == 3 ? nullptr : std::make_unique(), + i == 4 ? nullptr : std::make_shared(), + this->angular_flux_storage_map_); + }); + } +} + +TYPED_TEST(FormulationUpdaterDriftDiffusionTest, UpdateFixedTermTest) { + system::EnergyGroup group_number(this->group_number); + quadrature::QuadraturePointIndex angle_index(this->angle_index); + bart::system::Index scalar_index{this->group_number, 0}; + + EXPECT_CALL(*this->mock_lhs_obs_ptr_, GetFixedTermPtr(scalar_index)).WillOnce(DoDefault()); + EXPECT_CALL(*this->mock_rhs_obs_ptr_, GetFixedTermPtr(scalar_index)).WillOnce(DoDefault()); + EXPECT_CALL(*this->integrated_flux_calculator_obs_ptr_, Integrate(ContainerEq(this->group_angular_flux_))) + .WillOnce(Return(this->integrated_angular_flux_)); + std::array moment_index{this->group_number, 0, 0}; + EXPECT_CALL(*this->high_order_moments_ptr_, GetMoment(moment_index)).WillOnce(ReturnRef(this->group_scalar_flux_)); + + for (auto& cell : this->cells_) { + EXPECT_CALL(*this->diffusion_formulation_obs_ptr_, FillCellStreamingTerm(_, cell, this->group_number)); + EXPECT_CALL(*this->diffusion_formulation_obs_ptr_, FillCellCollisionTerm(_, cell, this->group_number)); + EXPECT_CALL(*this->diffusion_formulation_obs_ptr_, FillCellFixedSource(_, cell, this->group_number)); + EXPECT_CALL(*this->drift_diffusion_formulation_obs_ptr_, + FillCellDriftDiffusionTerm(_, + cell, + system::EnergyGroup(this->group_number), + this->group_scalar_flux_, + this->integrated_angular_flux_)); + if (cell->at_boundary()) { + int faces_per_cell = dealii::GeometryInfodim>::faces_per_cell; + for (int face = 0; face < faces_per_cell; ++face) { + if (cell->face(face)->at_boundary()) { + problem::Boundary boundary_id = static_cast(cell->face(face)->boundary_id()); + + using BoundaryType = typename formulation::scalar::DiffusionIdim>::BoundaryType; + BoundaryType boundary_type = BoundaryType::kVacuum; + if (this->reflective_boundaries.count(boundary_id) == 1) + boundary_type = BoundaryType::kReflective; + + EXPECT_CALL(*this->diffusion_formulation_obs_ptr_, FillBoundaryTerm(_, cell, face, boundary_type)); + } + } + } + } + + EXPECT_CALL(*this->stamper_obs_ptr_, StampMatrix(Ref(*this->matrix_to_stamp), _)) + .Times(3) + .WillRepeatedly(DoDefault()); + EXPECT_CALL(*this->stamper_obs_ptr_, StampBoundaryMatrix(Ref(*this->matrix_to_stamp), _)).WillOnce(DoDefault()); + EXPECT_CALL(*this->stamper_obs_ptr_, StampVector(Ref(*this->vector_to_stamp), _)).WillOnce(DoDefault()); + + this->test_updater_ptr_->UpdateFixedTerms(this->test_system_, group_number, angle_index); + EXPECT_TRUE(test_helpers::AreEqual(this->expected_result, *this->matrix_to_stamp)); + EXPECT_TRUE(test_helpers::AreEqual(this->expected_vector_result, *this->vector_to_stamp)); +} + +} // namespace diff --git a/src/framework/builder/framework_builder.cpp b/src/framework/builder/framework_builder.cpp index 163ae307f..85b7996e1 100644 --- a/src/framework/builder/framework_builder.cpp +++ b/src/framework/builder/framework_builder.cpp @@ -24,7 +24,7 @@ #include "formulation/scalar/diffusion.h" #include "formulation/stamper.h" #include "formulation/updater/saaf_updater.h" -#include "formulation/updater/diffusion_updater.h" +#include "formulation/updater/diffusion_updater.hpp" // Framework class #include "framework/framework.hpp" diff --git a/src/framework/builder/tests/framework_builder_integration_tests.cpp b/src/framework/builder/tests/framework_builder_integration_tests.cpp index b20211b98..a20df8395 100644 --- a/src/framework/builder/tests/framework_builder_integration_tests.cpp +++ b/src/framework/builder/tests/framework_builder_integration_tests.cpp @@ -18,7 +18,7 @@ #include "formulation/scalar/diffusion.h" #include "formulation/angular/self_adjoint_angular_flux.h" #include "formulation/updater/saaf_updater.h" -#include "formulation/updater/diffusion_updater.h" +#include "formulation/updater/diffusion_updater.hpp" #include "formulation/stamper.h" #include "instrumentation/instrument.h" #include "instrumentation/basic_instrument.h" diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp new file mode 100644 index 000000000..afbc4dbe6 --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp @@ -0,0 +1,35 @@ +#include "quadrature/calculators/drift_diffusion_integrated_flux.hpp" + +namespace bart::quadrature::calculators { + +template +DriftDiffusionIntegratedFlux::DriftDiffusionIntegratedFlux(std::shared_ptr quadrature_set_ptr) + : quadrature_set_ptr_(quadrature_set_ptr) { + this->AssertPointerNotNull(quadrature_set_ptr_.get(), "quadrature set", "DriftDiffusionIntegratedFlux constructor"); +} + +template +auto DriftDiffusionIntegratedFlux::Integrate(const VectorMap& vector_map) const -> Vector { + using Index = quadrature::QuadraturePointIndex; + const VectorMap::size_type n_quadrature_points{ this->quadrature_set_ptr_->size() }; + AssertThrow(vector_map.size() == n_quadrature_points, dealii::ExcMessage("Error in DriftDiffusionIntegratedFlux " + "integration, angular flux map is not the " + "same size as the quadrature set")); + const Vector::size_type vector_size{ vector_map.begin()->second->size()}; + Vector result_vector(vector_size); + + for (VectorMap::size_type i = 0; i < n_quadrature_points; ++i) { + auto& quadrature_point = *quadrature_set_ptr_->GetQuadraturePoint(Index(i)); + const double weight{ quadrature_point.weight() }; + const auto position{ quadrature_point.cartesian_position_tensor() }; + result_vector.add(weight * position * position, *vector_map.at(Index(i))); + } + + return result_vector; +} + +template class DriftDiffusionIntegratedFlux<1>; +template class DriftDiffusionIntegratedFlux<2>; +template class DriftDiffusionIntegratedFlux<3>; + +} // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp new file mode 100644 index 000000000..4559c9c19 --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp @@ -0,0 +1,28 @@ +#ifndef BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_HPP_ +#define BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_HPP_ + +#include "quadrature/calculators/drift_diffusion_integrated_flux_i.hpp" + +#include "quadrature/quadrature_set_i.h" +#include "utility/has_dependencies.h" + +namespace bart::quadrature::calculators { + +template +class DriftDiffusionIntegratedFlux : public DriftDiffusionIntegratedFluxI, public utility::HasDependencies { + public: + using QuadratureSet = typename quadrature::QuadratureSetI; + + DriftDiffusionIntegratedFlux(std::shared_ptr); + + [[nodiscard]] auto Integrate(const VectorMap&) const -> Vector override; + + auto quadrature_set_ptr() const -> QuadratureSet* { return quadrature_set_ptr_.get(); } + + private: + std::shared_ptr quadrature_set_ptr_{ nullptr }; +}; + +} // namespace bart::quadrature::calculators + +#endif //BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_HPP_ diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp new file mode 100644 index 000000000..5c20c7737 --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp @@ -0,0 +1,24 @@ +#ifndef BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_I_HPP_ +#define BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_I_HPP_ + +#include + +#include + +#include "quadrature/quadrature_types.h" + +namespace bart::quadrature::calculators { + +class DriftDiffusionIntegratedFluxI { + public: + virtual ~DriftDiffusionIntegratedFluxI() = default; + using Vector = dealii::Vector; + using VectorPtr = std::shared_ptr>; + using VectorMap = std::map; + + virtual auto Integrate(const VectorMap&) const -> Vector = 0; +}; + +} // namespace bart::quadrature::calculators + +#endif //BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_I_HPP_ diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp new file mode 100644 index 000000000..985c20a1e --- /dev/null +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp @@ -0,0 +1,19 @@ +#ifndef BART_SRC_QUADRATURE_CALCULATORS_TESTS_DRIFT_DIFFUSION_INTEGRATED_FLUX_MOCK_HPP_ +#define BART_SRC_QUADRATURE_CALCULATORS_TESTS_DRIFT_DIFFUSION_INTEGRATED_FLUX_MOCK_HPP_ + +#include "quadrature/calculators/drift_diffusion_integrated_flux_i.hpp" + +#include "test_helpers/gmock_wrapper.h" + +namespace bart::quadrature::calculators { + +class DriftDiffusionIntegratedFluxMock : public DriftDiffusionIntegratedFluxI { + public: + MOCK_METHOD(Vector, Integrate, (const VectorMap&), (const, override)); +}; + +} // namespace bart::quadrature::calculators + + + +#endif //BART_SRC_QUADRATURE_CALCULATORS_TESTS_DRIFT_DIFFUSION_INTEGRATED_FLUX_MOCK_HPP_ diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp new file mode 100644 index 000000000..275c2b5f3 --- /dev/null +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp @@ -0,0 +1,129 @@ +#include "quadrature/calculators/drift_diffusion_integrated_flux.hpp" + +#include + +#include "quadrature/tests/quadrature_set_mock.h" +#include "quadrature/tests/quadrature_point_mock.h" +#include "quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp" +#include "test_helpers/gmock_wrapper.h" +#include "test_helpers/test_assertions.hpp" + +namespace { + +using namespace bart; + +using ::testing::NiceMock, ::testing::Return, ::testing::DoDefault, ::testing::ContainerEq; + +template +class DriftDiffusionIntegratedFluxTest : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; + using QuadraturePointType = NiceMock>; + using QuadratureSetType = NiceMock>; + using TestIntegrator = typename quadrature::calculators::DriftDiffusionIntegratedFlux; + using Vector = dealii::Vector; + using VectorPtr = std::shared_ptr>; + using VectorMap = std::map; + + DriftDiffusionIntegratedFluxTest() + : expected_result_(expected_result_values_.cbegin(), expected_result_values_.cend()) {} + + // Test parameters + static constexpr int n_quadrature_points{ 3 }; + static constexpr int n_total_dofs{ 2 }; + + // Test object + std::unique_ptr test_integrator_{ nullptr }; + + // Mock objects + std::shared_ptr quadrature_set_ptr_{ nullptr }; + std::array, n_quadrature_points> mock_quadrature_points_; + + // Supporting objects + VectorMap angular_flux_map_{}; + const std::vector expected_result_values_{5501.5 * dim, 11003 * dim}; + const Vector expected_result_; + + auto SetUp() -> void override; +}; + +template +auto DriftDiffusionIntegratedFluxTest::SetUp() -> void { + quadrature_set_ptr_ = std::make_shared(); + test_integrator_ = std::move(std::make_unique(quadrature_set_ptr_)); + + using Index = quadrature::QuadraturePointIndex; + const std::array weights{0.15, 0.25, 0.6}; + for (int i = 0; i < n_quadrature_points; ++i) { + auto quadrature_point_ptr = std::make_shared(); + mock_quadrature_points_.at(i) = quadrature_point_ptr; + dealii::Tensor<1, dim> position_tensor; + for (int j = 0; j < dim; ++j) + position_tensor[j] = i + 1; + ON_CALL(*quadrature_point_ptr, weight()).WillByDefault(Return(weights.at(i))); + ON_CALL(*quadrature_point_ptr, cartesian_position_tensor()).WillByDefault(Return(position_tensor)); + ON_CALL(*quadrature_set_ptr_, GetQuadraturePoint(Index(i))).WillByDefault(Return(quadrature_point_ptr)); + auto vector_ptr = std::make_shared(n_total_dofs); + for (int j = 0; j < n_total_dofs; ++j) + (*vector_ptr)[j] = std::pow(10.0, i + 1) * (j + 1); + angular_flux_map_.insert({Index(i), vector_ptr}); + } + ON_CALL(*quadrature_set_ptr_, size()).WillByDefault(Return(n_quadrature_points)); +} + +TYPED_TEST_SUITE(DriftDiffusionIntegratedFluxTest, bart::testing::AllDimensions); + +TYPED_TEST(DriftDiffusionIntegratedFluxTest, Constructor) { + constexpr int dim = this->dim; + using TestIntegrator = typename quadrature::calculators::DriftDiffusionIntegratedFlux; + using QuadratureSetType = typename quadrature::QuadratureSetMock; + + auto quadrature_set_ptr = std::make_shared(); + TestIntegrator integrator(quadrature_set_ptr); + ASSERT_NE(nullptr, integrator.quadrature_set_ptr()); + EXPECT_EQ(integrator.quadrature_set_ptr(), quadrature_set_ptr.get()); +} + +TYPED_TEST(DriftDiffusionIntegratedFluxTest, ConstructorBadDependencies) { + constexpr int dim = this->dim; + using TestIntegrator = typename quadrature::calculators::DriftDiffusionIntegratedFlux; + + EXPECT_ANY_THROW({ + TestIntegrator integrator(nullptr); + }); +} + +TYPED_TEST(DriftDiffusionIntegratedFluxTest, Integrate) { + EXPECT_CALL(*this->quadrature_set_ptr_, size()).WillOnce(DoDefault()); + for (int i = 0; i < this->n_quadrature_points; ++i) { + using Index = quadrature::QuadraturePointIndex; + EXPECT_CALL(*this->quadrature_set_ptr_, GetQuadraturePoint(Index(i))).WillOnce(DoDefault()); + } + for (auto& quadrature_point : this->mock_quadrature_points_) { + EXPECT_CALL(*quadrature_point, weight()).WillOnce(DoDefault()); + EXPECT_CALL(*quadrature_point, cartesian_position_tensor()).WillOnce(DoDefault()); + } + auto result = this->test_integrator_->Integrate(this->angular_flux_map_); + EXPECT_TRUE(test_helpers::AreEqual(result, this->expected_result_)); +} + +TYPED_TEST(DriftDiffusionIntegratedFluxTest, IntegrateBadAngularFluxSize) { + using Vector = dealii::Vector; + using VectorPtr = std::shared_ptr>; + using VectorMap = std::map; + using Index = quadrature::QuadraturePointIndex; + + EXPECT_CALL(*this->quadrature_set_ptr_, size()).WillOnce(DoDefault()); + + VectorMap bad_angular_flux_map; + bad_angular_flux_map.insert({Index(0), std::make_shared()}); + EXPECT_ANY_THROW({ + [[maybe_unused]] auto result = this->test_integrator_->Integrate(bad_angular_flux_map); + }); +} + + + + + +} // namespace diff --git a/src/test_helpers/test_assertions.cpp b/src/test_helpers/test_assertions.cpp index cf3886be5..1cd1ad0c0 100644 --- a/src/test_helpers/test_assertions.cpp +++ b/src/test_helpers/test_assertions.cpp @@ -1,5 +1,6 @@ #include "test_helpers/test_assertions.hpp" +#include #include namespace bart { @@ -12,6 +13,9 @@ using ::testing::AssertionSuccess; using DealiiVector = dealii::Vector; using FullMatrix = dealii::FullMatrix; + +template +using Tensor = dealii::Tensor<1, dim>; } template <> @@ -33,6 +37,33 @@ auto AreEqual>(const DealiiVector& expected, const Dealii return AssertionSuccess(); } +template <> +auto AreEqual>(const Tensor<1>& expected, const Tensor<1>& result, const double tol) -> AssertionResult { + if (abs(result[0] - expected[0]) > tol) + return AssertionFailure() << "Entry (" << 0 << ") has value: " << result[0] << ", expected: " << expected[0]; + return AssertionSuccess(); +} + +template <> +auto AreEqual>(const Tensor<2>& expected, const Tensor<2>& result, const double tol) -> AssertionResult { + for (int i = 0; i < 2; ++i) { + if (abs(result[i] - expected[i]) > tol) { + return AssertionFailure() << "Entry (" << i << ") has value: " << result[i] << ", expected: " << expected[i]; + } + } + return AssertionSuccess(); +} + +template <> +auto AreEqual>(const Tensor<3>& expected, const Tensor<3>& result, const double tol) -> AssertionResult { + for (int i = 0; i < 3; ++i) { + if (abs(result[i] - expected[i]) > tol) { + return AssertionFailure() << "Entry (" << i << ") has value: " << result[i] << ", expected: " << expected[i]; + } + } + return AssertionSuccess(); +} + template <> auto AreEqual>(const std::vector& expected, const std::vector& result, const double tol) -> AssertionResult{ diff --git a/src/test_helpers/tests/test_assertions_test.cpp b/src/test_helpers/tests/test_assertions_test.cpp index 8056e9253..8a6005702 100644 --- a/src/test_helpers/tests/test_assertions_test.cpp +++ b/src/test_helpers/tests/test_assertions_test.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "test_helpers/dealii_test_domain.h" #include "test_helpers/test_helper_functions.h" @@ -123,6 +124,42 @@ TEST_F(TestAssertionsMatrixTests, GoodComparisonWithinTolerance) { EXPECT_TRUE(test_helpers::AreEqual(matrix_1, matrix_3, 1e-4)); } +template +class TestAssertionsTensorsAreEqual : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; + dealii::Tensor<1, dim> tensor_1_, tensor_2_; + + auto SetUp() -> void override; +}; + +template +auto TestAssertionsTensorsAreEqual::SetUp() -> void { + for (int i = 0; i < dim; ++i) { + tensor_1_[i] = test_helpers::RandomDouble(-100, 100); + tensor_2_[i] = test_helpers::RandomDouble(-100, 100); + } +} + +TYPED_TEST_SUITE(TestAssertionsTensorsAreEqual, bart::testing::AllDimensions); + +TYPED_TEST(TestAssertionsTensorsAreEqual, GoodComparison) { + EXPECT_TRUE(test_helpers::AreEqual(this->tensor_1_, this->tensor_1_)); + EXPECT_TRUE(test_helpers::AreEqual(this->tensor_2_, this->tensor_2_)); +} + +TYPED_TEST(TestAssertionsTensorsAreEqual, BadComparison) { + EXPECT_FALSE(test_helpers::AreEqual(this->tensor_1_, this->tensor_2_)); + EXPECT_FALSE(test_helpers::AreEqual(this->tensor_2_, this->tensor_1_)); +} + +TYPED_TEST(TestAssertionsTensorsAreEqual, Tolerance) { + auto tensor_3 = this->tensor_1_; + tensor_3 *= (1 + 1e-8); + EXPECT_TRUE(test_helpers::AreEqual(this->tensor_1_, tensor_3)); + tensor_3 *= (1 + 1e-5); + EXPECT_FALSE(test_helpers::AreEqual(this->tensor_1_, tensor_3)); +} class TestAssertionsMPIMatricesTests : public ::testing::Test, public bart::testing::DealiiTestDomain<2> {