From 0e128fbb8e02c9ff24fafab042b7126d250252b9 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 7 Dec 2020 15:13:09 -0800 Subject: [PATCH 01/23] Initial commit of DriftDiffusionVectorCalculator. --- .../drift_diffusion_vector_calculator.hpp | 6 +++ .../drift_diffusion_vector_calculator_i.hpp | 40 +++++++++++++++++++ ...drift_diffusion_vector_calculator_test.cpp | 21 ++++++++++ 3 files changed, 67 insertions(+) create mode 100644 src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp create mode 100644 src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp create mode 100644 src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp 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..ed293ed2d --- /dev/null +++ b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp @@ -0,0 +1,6 @@ +#ifndef BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_HPP_ +#define BART_SRC_CALCULATOR_DRIFT_DIFFUSION_DRIFT_DIFFUSION_VECTOR_CALCULATOR_HPP_ + + + +#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..d0d49fe3c --- /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& shapte_gradient, + const double sigma_t, + const double diffusion_coefficient) = 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_test.cpp b/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp new file mode 100644 index 000000000..f26eb192d --- /dev/null +++ b/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp @@ -0,0 +1,21 @@ +#include "test_helpers/gmock_wrapper.h" + +#include "calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp" + +namespace { + +using namespace bart; + +template +class DriftDiffusionVectorCalculatorTest : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; +}; + +TYPED_TEST_SUITE(DriftDiffusionVectorCalculatorTest, bart::testing::AllDimensions); + +TYPED_TEST(DriftDiffusionVectorCalculatorTest, Dummy) { + EXPECT_TRUE(false); +} + +} // namespace From 8a6b5dd5581235cf280394db5903256cdeb2e6e5 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 7 Dec 2020 15:53:20 -0800 Subject: [PATCH 02/23] Added AreEqual specialization for dealii Tensors. --- src/test_helpers/test_assertions.cpp | 31 ++++++++++++++++ .../tests/test_assertions_test.cpp | 37 +++++++++++++++++++ 2 files changed, 68 insertions(+) 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> { From 7fb0bf4f568b3c90a81cf0740e233311d16a11ba Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 7 Dec 2020 16:08:36 -0800 Subject: [PATCH 03/23] Added implementation of DriftDiffusion calculation. --- .../drift_diffusion_vector_calculator.cpp | 28 ++++++ .../drift_diffusion_vector_calculator.hpp | 16 ++++ .../drift_diffusion_vector_calculator_i.hpp | 4 +- ...drift_diffusion_vector_calculator_test.cpp | 87 ++++++++++++++++++- 4 files changed, 131 insertions(+), 4 deletions(-) create mode 100644 src/calculator/drift_diffusion/drift_diffusion_vector_calculator.cpp 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 index ed293ed2d..55adc1daf 100644 --- a/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp +++ b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator.hpp @@ -1,6 +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 index d0d49fe3c..bfe351241 100644 --- a/src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp +++ b/src/calculator/drift_diffusion/drift_diffusion_vector_calculator_i.hpp @@ -30,9 +30,9 @@ class DriftDiffusionVectorCalculatorI { using Tensor = typename dealii::Tensor<1, dim>; virtual auto DriftDiffusion(const double scalar_flux, const double integrated_angular_flux, - const Tensor& shapte_gradient, + const Tensor& shape_gradient, const double sigma_t, - const double diffusion_coefficient) = 0; + const double diffusion_coefficient) const -> Tensor = 0; }; } // namespace bart::calculator::drift_diffusion 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 index f26eb192d..2984be5ad 100644 --- a/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp +++ b/src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_test.cpp @@ -1,21 +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, Dummy) { - EXPECT_TRUE(false); +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 From d100fd5d6cabd0fae0bc32efdf4f0bcdec516b8b Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 8 Dec 2020 09:48:54 -0800 Subject: [PATCH 04/23] Initial commit of DriftDiffusionIntegratedFlux classes. --- .../drift_diffusion_integrated_flux.hpp | 15 +++++++++++++ .../drift_diffusion_integrated_flux_i.hpp | 12 ++++++++++ .../drift_diffusion_integrated_flux_mock.hpp | 19 ++++++++++++++++ .../drift_diffusion_integrated_flux_test.cpp | 22 +++++++++++++++++++ 4 files changed, 68 insertions(+) create mode 100644 src/quadrature/calculators/drift_diffusion_integrated_flux.hpp create mode 100644 src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp create mode 100644 src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp create mode 100644 src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp 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..3ab5cce4c --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp @@ -0,0 +1,15 @@ +#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" + +namespace bart::quadrature::calculators { + +template +class DriftDiffusionIntegratedFlux : public DriftDiffusionIntegratedFluxI { + public: +}; + +} // 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..c3800fbc2 --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp @@ -0,0 +1,12 @@ +#ifndef BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_I_HPP_ +#define BART_SRC_QUADRATURE_CALCULATORS_DRIFT_DIFFUSION_INTEGRATED_FLUX_I_HPP_ + +namespace bart::quadrature::calculators { + +class DriftDiffusionIntegratedFluxI { + public: +}; + +} // 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..b86d898c0 --- /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: + +}; + +} // 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..aa2642d46 --- /dev/null +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp @@ -0,0 +1,22 @@ +#include "quadrature/calculators/drift_diffusion_integrated_flux.hpp" + +#include "quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp" +#include "test_helpers/gmock_wrapper.h" + +namespace { + +using namespace bart; + +template +class DriftDiffusionIntegratedFluxTest : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; +}; + +TYPED_TEST_SUITE(DriftDiffusionIntegratedFluxTest, bart::testing::AllDimensions); + +TYPED_TEST(DriftDiffusionIntegratedFluxTest, Dummy) { + EXPECT_TRUE(false); +} + +} // namespace From dc5774486b9e9f536b99cc317f4c8f1816c8e411 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 8 Dec 2020 10:02:45 -0800 Subject: [PATCH 05/23] Added constructor with quadrature set dependency. --- .../drift_diffusion_integrated_flux.cpp | 15 ++++++++ .../drift_diffusion_integrated_flux.hpp | 13 ++++++- .../drift_diffusion_integrated_flux_i.hpp | 1 + .../drift_diffusion_integrated_flux_test.cpp | 36 +++++++++++++++++-- 4 files changed, 62 insertions(+), 3 deletions(-) create mode 100644 src/quadrature/calculators/drift_diffusion_integrated_flux.cpp 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..26bcf7aac --- /dev/null +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp @@ -0,0 +1,15 @@ +#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 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 index 3ab5cce4c..d211e3876 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp @@ -3,11 +3,22 @@ #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 { +class DriftDiffusionIntegratedFlux : public DriftDiffusionIntegratedFluxI, public utility::HasDependencies { public: + using QuadratureSet = typename quadrature::QuadratureSetI; + + DriftDiffusionIntegratedFlux(std::shared_ptr); + + auto quadrature_set_ptr() const -> QuadratureSet* { return quadrature_set_ptr_.get(); } + + private: + std::shared_ptr quadrature_set_ptr_{ nullptr }; }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp index c3800fbc2..c0e0fce45 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp @@ -5,6 +5,7 @@ namespace bart::quadrature::calculators { class DriftDiffusionIntegratedFluxI { public: + virtual ~DriftDiffusionIntegratedFluxI() = default; }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp index aa2642d46..8307ea14b 100644 --- a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp @@ -1,5 +1,8 @@ #include "quadrature/calculators/drift_diffusion_integrated_flux.hpp" +#include + +#include "quadrature/tests/quadrature_set_mock.h" #include "quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp" #include "test_helpers/gmock_wrapper.h" @@ -11,12 +14,41 @@ template class DriftDiffusionIntegratedFluxTest : public ::testing::Test { public: static constexpr int dim = DimensionWrapper::value; + using QuadratureSetType = typename quadrature::QuadratureSetMock; + + // Mock objects + std::shared_ptr quadrature_set_ptr_; + + auto SetUp() -> void override; }; +template +auto DriftDiffusionIntegratedFluxTest::SetUp() -> void { + quadrature_set_ptr_ = std::make_shared(); +} + TYPED_TEST_SUITE(DriftDiffusionIntegratedFluxTest, bart::testing::AllDimensions); -TYPED_TEST(DriftDiffusionIntegratedFluxTest, Dummy) { - EXPECT_TRUE(false); +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); + }); } + + } // namespace From 33f202e9c83eba051263e6c66ed82c3fcd7adf60 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 8 Dec 2020 11:33:18 -0800 Subject: [PATCH 06/23] Added implementation and testing of Integrate. --- .../drift_diffusion_integrated_flux.cpp | 24 ++++++ .../drift_diffusion_integrated_flux.hpp | 2 + .../drift_diffusion_integrated_flux_i.hpp | 11 +++ .../drift_diffusion_integrated_flux_mock.hpp | 2 +- .../drift_diffusion_integrated_flux_test.cpp | 75 ++++++++++++++++++- 5 files changed, 111 insertions(+), 3 deletions(-) diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp index 26bcf7aac..6f53a9071 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp @@ -8,6 +8,30 @@ DriftDiffusionIntegratedFlux::DriftDiffusionIntegratedFlux(std::shared_ptr< this->AssertPointerNotNull(quadrature_set_ptr_.get(), "quadrature set", "DriftDiffusionIntegratedFlux constructor"); } +template +auto DriftDiffusionIntegratedFlux::Integrate(const VectorMap& vector_map) const -> std::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))); + } + + std::vector return_vector(vector_size); + for (Vector::size_type i = 0; i < vector_size; ++i) + return_vector.at(i) = result_vector[i]; + + return return_vector; +} + template class DriftDiffusionIntegratedFlux<1>; template class DriftDiffusionIntegratedFlux<2>; template class DriftDiffusionIntegratedFlux<3>; diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp index d211e3876..9112e7f92 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp @@ -15,6 +15,8 @@ class DriftDiffusionIntegratedFlux : public DriftDiffusionIntegratedFluxI, publi DriftDiffusionIntegratedFlux(std::shared_ptr); + [[nodiscard]] auto Integrate(const VectorMap&) const -> std::vector override; + auto quadrature_set_ptr() const -> QuadratureSet* { return quadrature_set_ptr_.get(); } private: diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp index c0e0fce45..6dbf86397 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp @@ -1,11 +1,22 @@ #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 -> std::vector = 0; }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp index b86d898c0..31b60521c 100644 --- a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp @@ -9,7 +9,7 @@ namespace bart::quadrature::calculators { class DriftDiffusionIntegratedFluxMock : public DriftDiffusionIntegratedFluxI { public: - + MOCK_METHOD(std::vector, Integrate, (const VectorMap&), (const, override)); }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp index 8307ea14b..aef372d5b 100644 --- a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp @@ -3,21 +3,42 @@ #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 QuadratureSetType = typename quadrature::QuadratureSetMock; + 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; + + // 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_; + std::shared_ptr quadrature_set_ptr_{ nullptr }; + std::array, n_quadrature_points> mock_quadrature_points_; + + // Supporting objects + VectorMap angular_flux_map_{}; + std::vector expected_result_{5501.5 * dim, 11003 * dim}; auto SetUp() -> void override; }; @@ -25,6 +46,25 @@ class DriftDiffusionIntegratedFluxTest : public ::testing::Test { 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); @@ -49,6 +89,37 @@ TYPED_TEST(DriftDiffusionIntegratedFluxTest, ConstructorBadDependencies) { }); } +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 From fbf12ead8ffd0821e74426a54ccf2f108434c4d1 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 8 Dec 2020 14:28:08 -0800 Subject: [PATCH 07/23] Updated Integrate to return a dealii::Vector for compatibility to ValueAtQuadrature. --- .../calculators/drift_diffusion_integrated_flux.cpp | 8 ++------ .../calculators/drift_diffusion_integrated_flux.hpp | 2 +- .../calculators/drift_diffusion_integrated_flux_i.hpp | 2 +- .../tests/drift_diffusion_integrated_flux_mock.hpp | 2 +- .../tests/drift_diffusion_integrated_flux_test.cpp | 6 +++++- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp index 6f53a9071..afbc4dbe6 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.cpp @@ -9,7 +9,7 @@ DriftDiffusionIntegratedFlux::DriftDiffusionIntegratedFlux(std::shared_ptr< } template -auto DriftDiffusionIntegratedFlux::Integrate(const VectorMap& vector_map) const -> std::vector { +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 " @@ -25,11 +25,7 @@ auto DriftDiffusionIntegratedFlux::Integrate(const VectorMap& vector_map) c result_vector.add(weight * position * position, *vector_map.at(Index(i))); } - std::vector return_vector(vector_size); - for (Vector::size_type i = 0; i < vector_size; ++i) - return_vector.at(i) = result_vector[i]; - - return return_vector; + return result_vector; } template class DriftDiffusionIntegratedFlux<1>; diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp index 9112e7f92..4559c9c19 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux.hpp @@ -15,7 +15,7 @@ class DriftDiffusionIntegratedFlux : public DriftDiffusionIntegratedFluxI, publi DriftDiffusionIntegratedFlux(std::shared_ptr); - [[nodiscard]] auto Integrate(const VectorMap&) const -> std::vector override; + [[nodiscard]] auto Integrate(const VectorMap&) const -> Vector override; auto quadrature_set_ptr() const -> QuadratureSet* { return quadrature_set_ptr_.get(); } diff --git a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp index 6dbf86397..5c20c7737 100644 --- a/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp +++ b/src/quadrature/calculators/drift_diffusion_integrated_flux_i.hpp @@ -16,7 +16,7 @@ class DriftDiffusionIntegratedFluxI { using VectorPtr = std::shared_ptr>; using VectorMap = std::map; - virtual auto Integrate(const VectorMap&) const -> std::vector = 0; + virtual auto Integrate(const VectorMap&) const -> Vector = 0; }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp index 31b60521c..985c20a1e 100644 --- a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_mock.hpp @@ -9,7 +9,7 @@ namespace bart::quadrature::calculators { class DriftDiffusionIntegratedFluxMock : public DriftDiffusionIntegratedFluxI { public: - MOCK_METHOD(std::vector, Integrate, (const VectorMap&), (const, override)); + MOCK_METHOD(Vector, Integrate, (const VectorMap&), (const, override)); }; } // namespace bart::quadrature::calculators diff --git a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp index aef372d5b..275c2b5f3 100644 --- a/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp +++ b/src/quadrature/calculators/tests/drift_diffusion_integrated_flux_test.cpp @@ -25,6 +25,9 @@ class DriftDiffusionIntegratedFluxTest : public ::testing::Test { 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 }; @@ -38,7 +41,8 @@ class DriftDiffusionIntegratedFluxTest : public ::testing::Test { // Supporting objects VectorMap angular_flux_map_{}; - std::vector expected_result_{5501.5 * dim, 11003 * dim}; + const std::vector expected_result_values_{5501.5 * dim, 11003 * dim}; + const Vector expected_result_; auto SetUp() -> void override; }; From 804cff024a7f56231fe414b4f76aaa994f3dd99b Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Wed, 9 Dec 2020 08:41:14 -0800 Subject: [PATCH 08/23] Initial commit of formulation::scalar::DriftDiffusion --- src/formulation/scalar/drift_diffusion.cpp | 9 ++++++++ src/formulation/scalar/drift_diffusion.hpp | 15 +++++++++++++ src/formulation/scalar/drift_diffusion_i.hpp | 14 +++++++++++++ .../scalar/tests/drift_diffusion_test.cpp | 21 +++++++++++++++++++ 4 files changed, 59 insertions(+) create mode 100644 src/formulation/scalar/drift_diffusion.cpp create mode 100644 src/formulation/scalar/drift_diffusion.hpp create mode 100644 src/formulation/scalar/drift_diffusion_i.hpp create mode 100644 src/formulation/scalar/tests/drift_diffusion_test.cpp diff --git a/src/formulation/scalar/drift_diffusion.cpp b/src/formulation/scalar/drift_diffusion.cpp new file mode 100644 index 000000000..2037d8e78 --- /dev/null +++ b/src/formulation/scalar/drift_diffusion.cpp @@ -0,0 +1,9 @@ +#include "formulation/scalar/drift_diffusion.hpp" + +namespace bart::formulation::scalar { + +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..f39d26469 --- /dev/null +++ b/src/formulation/scalar/drift_diffusion.hpp @@ -0,0 +1,15 @@ +#ifndef BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_HPP_ +#define BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_HPP_ + +#include "formulation/scalar/drift_diffusion_i.hpp" + +namespace bart::formulation::scalar { + +template +class DriftDiffusion : public DriftDiffusionI { + public: +}; + +} // 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..20e7f8fd0 --- /dev/null +++ b/src/formulation/scalar/drift_diffusion_i.hpp @@ -0,0 +1,14 @@ +#ifndef BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_HPP_ +#define BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_HPP_ + +namespace bart::formulation::scalar { + +template +class DriftDiffusionI { + public: + virtual ~DriftDiffusionI() = default; +}; + +} // namespace bart::formulation::scalar + +#endif //BART_SRC_FORMULATION_SCALAR_DRIFT_DIFFUSION_I_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..a4acc12a2 --- /dev/null +++ b/src/formulation/scalar/tests/drift_diffusion_test.cpp @@ -0,0 +1,21 @@ +#include "formulation/scalar/drift_diffusion.hpp" + +#include "test_helpers/gmock_wrapper.h" + +namespace { + +using namespace bart; + +template +class DriftDiffusionFormulationTest : public ::testing::Test { + public: + static constexpr int dim = DimensionWrapper::value; +}; + +TYPED_TEST_SUITE(DriftDiffusionFormulationTest, bart::testing::AllDimensions); + +TYPED_TEST(DriftDiffusionFormulationTest, Dummy) { + EXPECT_TRUE(false); +} + +} // namespace From 7094f59590980d88cb4c63a57798b3ffe4c34555 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 09:08:47 -0800 Subject: [PATCH 09/23] Added mock DriftDiffusionCalculator. --- ...drift_diffusion_vector_calculator_mock.hpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 src/calculator/drift_diffusion/tests/drift_diffusion_vector_calculator_mock.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_ From d06a800b2a7f47428357971ef55eb79ccd6fd35c Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 09:40:55 -0800 Subject: [PATCH 10/23] Added constructor and dependencies. --- src/formulation/scalar/drift_diffusion.cpp | 13 ++ src/formulation/scalar/drift_diffusion.hpp | 25 +++- .../scalar/tests/drift_diffusion_test.cpp | 136 +++++++++++++++++- 3 files changed, 171 insertions(+), 3 deletions(-) diff --git a/src/formulation/scalar/drift_diffusion.cpp b/src/formulation/scalar/drift_diffusion.cpp index 2037d8e78..bbe650e28 100644 --- a/src/formulation/scalar/drift_diffusion.cpp +++ b/src/formulation/scalar/drift_diffusion.cpp @@ -2,6 +2,19 @@ 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); +} + template class DriftDiffusion<1>; template class DriftDiffusion<2>; template class DriftDiffusion<3>; diff --git a/src/formulation/scalar/drift_diffusion.hpp b/src/formulation/scalar/drift_diffusion.hpp index f39d26469..b36bfb89e 100644 --- a/src/formulation/scalar/drift_diffusion.hpp +++ b/src/formulation/scalar/drift_diffusion.hpp @@ -3,11 +3,34 @@ #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 { +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; + + DriftDiffusion(std::shared_ptr, + std::shared_ptr, + std::shared_ptr); + + 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 }; }; } // namespace bart::formulation::scalar diff --git a/src/formulation/scalar/tests/drift_diffusion_test.cpp b/src/formulation/scalar/tests/drift_diffusion_test.cpp index a4acc12a2..c663913a5 100644 --- a/src/formulation/scalar/tests/drift_diffusion_test.cpp +++ b/src/formulation/scalar/tests/drift_diffusion_test.cpp @@ -1,21 +1,153 @@ #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" namespace { using namespace bart; +using test_helpers::AreEqual; + +using ::testing::NiceMock, ::testing::Return, ::testing::DoDefault, ::testing::_; + 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 FiniteElement = NiceMock>; + using Material = NiceMock; + + DriftDiffusionFormulationTest() : dof_handler_(triangulation_), finite_element_(1) {}; + + 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 }; + dealii::FullMatrix expected_result_q_0_, expected_result_q_1_; + + auto SetUp() -> void override; }; +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_)); + + auto ReturnShapeGradient = [&](const int i, const int q) { + 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; + }; + + 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) { + ON_CALL(*finite_element_mock_ptr_, ShapeValue(i,q)).WillByDefault(Return(1 + i + q)); + ON_CALL(*finite_element_mock_ptr_, ShapeGradient(i,q)).WillByDefault(ReturnShapeGradient); + } + } + + 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)); + + std::unordered_map> sigma_t{{material_id_, {1.0, 2.0}}}; + std::unordered_map> diffusion_coef{{material_id_, {0.5, 1.0}}}; + ON_CALL(mock_material, GetSigT()).WillByDefault(Return(sigma_t)); + ON_CALL(mock_material, GetDiffusionCoef()).WillByDefault(Return(sigma_t)); + 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}; + expected_result_q_0_ = dealii::FullMatrix(2, 2, expected_result_q_0_values.begin()); + expected_result_q_1_ = dealii::FullMatrix(2, 2, expected_result_q_1_values.begin()); +} + TYPED_TEST_SUITE(DriftDiffusionFormulationTest, bart::testing::AllDimensions); -TYPED_TEST(DriftDiffusionFormulationTest, Dummy) { - EXPECT_TRUE(false); +TYPED_TEST(DriftDiffusionFormulationTest, ConstructorAndDependencyGetters) { + 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 < 3; ++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); + }); + } } + + } // namespace From a659f1c5e2868759ead1a966aee4d65065d47894 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 10:02:55 -0800 Subject: [PATCH 11/23] Updated call to ValueAtQuadrature to use a reference. --- src/domain/finite_element/finite_element.cc | 3 +-- src/domain/finite_element/finite_element.h | 3 +-- src/domain/finite_element/finite_element_i.h | 3 +-- src/domain/finite_element/tests/finite_element_mock.h | 2 +- 4 files changed, 4 insertions(+), 7 deletions(-) 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)); From 1c0951ea2bc7950cbd3982a235aea3bd1f3f7f60 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 10:57:12 -0800 Subject: [PATCH 12/23] Added FillCellDriftDiffusionTerm. --- src/formulation/scalar/drift_diffusion.cpp | 32 ++++ src/formulation/scalar/drift_diffusion.hpp | 10 ++ src/formulation/scalar/drift_diffusion_i.hpp | 13 ++ .../scalar/tests/drift_diffusion_test.cpp | 155 ++++++++++++++---- 4 files changed, 177 insertions(+), 33 deletions(-) diff --git a/src/formulation/scalar/drift_diffusion.cpp b/src/formulation/scalar/drift_diffusion.cpp index bbe650e28..b9f8af9a8 100644 --- a/src/formulation/scalar/drift_diffusion.cpp +++ b/src/formulation/scalar/drift_diffusion.cpp @@ -13,6 +13,38 @@ DriftDiffusion::DriftDiffusion(std::shared_ptr finite_elemen 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 { + finite_element_ptr_->SetCell(cell_ptr); + const int 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>; diff --git a/src/formulation/scalar/drift_diffusion.hpp b/src/formulation/scalar/drift_diffusion.hpp index b36bfb89e..92a5221a6 100644 --- a/src/formulation/scalar/drift_diffusion.hpp +++ b/src/formulation/scalar/drift_diffusion.hpp @@ -18,10 +18,18 @@ class DriftDiffusion : public DriftDiffusionI, public utility::HasDependenc 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(); } @@ -31,6 +39,8 @@ class DriftDiffusion : public DriftDiffusionI, public utility::HasDependenc 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 diff --git a/src/formulation/scalar/drift_diffusion_i.hpp b/src/formulation/scalar/drift_diffusion_i.hpp index 20e7f8fd0..b4f6c562f 100644 --- a/src/formulation/scalar/drift_diffusion_i.hpp +++ b/src/formulation/scalar/drift_diffusion_i.hpp @@ -1,12 +1,25 @@ #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 diff --git a/src/formulation/scalar/tests/drift_diffusion_test.cpp b/src/formulation/scalar/tests/drift_diffusion_test.cpp index c663913a5..2dacbdef3 100644 --- a/src/formulation/scalar/tests/drift_diffusion_test.cpp +++ b/src/formulation/scalar/tests/drift_diffusion_test.cpp @@ -11,6 +11,7 @@ #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 { @@ -18,7 +19,7 @@ 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 { @@ -27,11 +28,15 @@ class DriftDiffusionFormulationTest : public ::testing::Test { 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 }; @@ -46,11 +51,33 @@ class DriftDiffusionFormulationTest : public ::testing::Test { const int dofs_per_cell_{ 2 }; const int cell_quadrature_points_{ 2 }; const int material_id_{ 1 }; - dealii::FullMatrix expected_result_q_0_, expected_result_q_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(); @@ -60,26 +87,12 @@ auto DriftDiffusionFormulationTest::SetUp() -> void { 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_)); - auto ReturnShapeGradient = [&](const int i, const int q) { - 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; - }; - 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(ReturnShapeGradient); + ON_CALL(*finite_element_mock_ptr_, ShapeGradient(i,q)).WillByDefault(Return(shape_gradient)); } } @@ -96,10 +109,8 @@ auto DriftDiffusionFormulationTest::SetUp() -> void { ON_CALL(*drift_diffusion_calculator_mock_ptr_, DriftDiffusion(_, _, _, _, _)).WillByDefault(Return(drift_diffusion)); - std::unordered_map> sigma_t{{material_id_, {1.0, 2.0}}}; - std::unordered_map> diffusion_coef{{material_id_, {0.5, 1.0}}}; - ON_CALL(mock_material, GetSigT()).WillByDefault(Return(sigma_t)); - ON_CALL(mock_material, GetDiffusionCoef()).WillByDefault(Return(sigma_t)); + 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); @@ -112,21 +123,42 @@ auto DriftDiffusionFormulationTest::SetUp() -> void { } } - 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}; - expected_result_q_0_ = dealii::FullMatrix(2, 2, expected_result_q_0_values.begin()); - expected_result_q_1_ = dealii::FullMatrix(2, 2, expected_result_q_1_values.begin()); +// 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_, @@ -138,7 +170,7 @@ TYPED_TEST(DriftDiffusionFormulationTest, ConstructorAndDependencyGetters) { TYPED_TEST(DriftDiffusionFormulationTest, ConstructorBadDependencies) { const int n_dependencies{ 3 }; - for (int i = 0; i < 3; ++i) { + 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, @@ -148,6 +180,63 @@ TYPED_TEST(DriftDiffusionFormulationTest, ConstructorBadDependencies) { } } +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 From 2f87165ff8e5e315f66a3ff0a6c43483e54bc6bb Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 13:45:41 -0800 Subject: [PATCH 13/23] Added Asserts to verify size of matrix to fill in DriftDiffusion. --- src/formulation/scalar/drift_diffusion.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/formulation/scalar/drift_diffusion.cpp b/src/formulation/scalar/drift_diffusion.cpp index b9f8af9a8..d48aaa581 100644 --- a/src/formulation/scalar/drift_diffusion.cpp +++ b/src/formulation/scalar/drift_diffusion.cpp @@ -22,8 +22,11 @@ auto DriftDiffusion::FillCellDriftDiffusionTerm(Matrix& to_fill, 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 int material_id{ cell_ptr->material_id() }; + 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()) }; From 726358b8995d2521e397c18651271e6f4180bac4 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 14:45:47 -0800 Subject: [PATCH 14/23] Renamed diffusion_updater.h and .cc to .hpp and .cpp --- src/formulation/factory/formulation_factories.h | 2 +- .../updater/{diffusion_updater.cc => diffusion_updater.cpp} | 2 +- .../updater/{diffusion_updater.h => diffusion_updater.hpp} | 1 + .../{diffusion_updater_test.cc => diffusion_updater_test.cpp} | 2 +- src/framework/builder/framework_builder.cpp | 2 +- .../builder/tests/framework_builder_integration_tests.cpp | 2 +- 6 files changed, 6 insertions(+), 5 deletions(-) rename src/formulation/updater/{diffusion_updater.cc => diffusion_updater.cpp} (99%) rename src/formulation/updater/{diffusion_updater.h => diffusion_updater.hpp} (97%) rename src/formulation/updater/tests/{diffusion_updater_test.cc => diffusion_updater_test.cpp} (99%) 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/updater/diffusion_updater.cc b/src/formulation/updater/diffusion_updater.cpp similarity index 99% rename from src/formulation/updater/diffusion_updater.cc rename to src/formulation/updater/diffusion_updater.cpp index fc332895b..6e309afbe 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 { diff --git a/src/formulation/updater/diffusion_updater.h b/src/formulation/updater/diffusion_updater.hpp similarity index 97% rename from src/formulation/updater/diffusion_updater.h rename to src/formulation/updater/diffusion_updater.hpp index a59aa0bef..200879431 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" 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/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" From c016e7c71a0689edac5f5365222063e3b328eaaa Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Thu, 10 Dec 2020 14:48:23 -0800 Subject: [PATCH 15/23] Changed stamper dependency in DiffusionUpdater to a shared ptr. --- src/formulation/updater/diffusion_updater.cpp | 2 +- src/formulation/updater/diffusion_updater.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/formulation/updater/diffusion_updater.cpp b/src/formulation/updater/diffusion_updater.cpp index 6e309afbe..80b78c566 100644 --- a/src/formulation/updater/diffusion_updater.cpp +++ b/src/formulation/updater/diffusion_updater.cpp @@ -12,7 +12,7 @@ DiffusionUpdater::DiffusionUpdater( std::unique_ptr stamper_ptr, std::unordered_set reflective_boundaries) : formulation_ptr_(std::move(formulation_ptr)), - stamper_ptr_(std::move(stamper_ptr)), + stamper_ptr_(std::shared_ptr(std::move(stamper_ptr))), reflective_boundaries_(reflective_boundaries) { AssertThrow(formulation_ptr_ != nullptr, dealii::ExcMessage("Error in constructor of DiffusionUpdater, " diff --git a/src/formulation/updater/diffusion_updater.hpp b/src/formulation/updater/diffusion_updater.hpp index 200879431..bf1907f4b 100644 --- a/src/formulation/updater/diffusion_updater.hpp +++ b/src/formulation/updater/diffusion_updater.hpp @@ -61,7 +61,7 @@ class DiffusionUpdater StamperType* stamper_ptr() const { return stamper_ptr_.get(); } private: std::unique_ptr formulation_ptr_; - std::unique_ptr stamper_ptr_; + std::shared_ptr stamper_ptr_; std::unordered_set reflective_boundaries_; }; From 9a67f5998253ea7aa971b195bea28d2e9a514ef1 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Fri, 11 Dec 2020 12:08:55 -0800 Subject: [PATCH 16/23] Added FixedUpdater class. --- src/formulation/updater/diffusion_updater.cpp | 62 ++++++++----------- src/formulation/updater/diffusion_updater.hpp | 17 ++--- src/formulation/updater/fixed_updater.hpp | 56 +++++++++++++++++ 3 files changed, 91 insertions(+), 44 deletions(-) create mode 100644 src/formulation/updater/fixed_updater.hpp diff --git a/src/formulation/updater/diffusion_updater.cpp b/src/formulation/updater/diffusion_updater.cpp index 80b78c566..6ccfd91e2 100644 --- a/src/formulation/updater/diffusion_updater.cpp +++ b/src/formulation/updater/diffusion_updater.cpp @@ -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::shared_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,35 @@ 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.hpp b/src/formulation/updater/diffusion_updater.hpp index bf1907f4b..415eb2952 100644 --- a/src/formulation/updater/diffusion_updater.hpp +++ b/src/formulation/updater/diffusion_updater.hpp @@ -23,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, @@ -59,6 +60,8 @@ 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::shared_ptr stamper_ptr_; diff --git a/src/formulation/updater/fixed_updater.hpp b/src/formulation/updater/fixed_updater.hpp new file mode 100644 index 000000000..aead90bf3 --- /dev/null +++ b/src/formulation/updater/fixed_updater.hpp @@ -0,0 +1,56 @@ +#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); + for (auto& vector_boundary_function : fixed_vector_boundary_functions_) + stamper_ptr_->StampBoundaryVector(*fixed_vector_ptr, vector_boundary_function); + } + + 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_ From 7c343e81b1185e6218b91cbee057aeb364b20f14 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 14 Dec 2020 09:51:10 -0800 Subject: [PATCH 17/23] Initial commit of DriftDiffusionUpdater. --- .../updater/drift_diffusion_updater.cpp | 11 ++++++++ .../updater/drift_diffusion_updater.hpp | 25 +++++++++++++++++++ .../tests/drift_diffusion_updater_test.cpp | 25 +++++++++++++++++++ 3 files changed, 61 insertions(+) create mode 100644 src/formulation/updater/drift_diffusion_updater.cpp create mode 100644 src/formulation/updater/drift_diffusion_updater.hpp create mode 100644 src/formulation/updater/tests/drift_diffusion_updater_test.cpp diff --git a/src/formulation/updater/drift_diffusion_updater.cpp b/src/formulation/updater/drift_diffusion_updater.cpp new file mode 100644 index 000000000..7ce72789f --- /dev/null +++ b/src/formulation/updater/drift_diffusion_updater.cpp @@ -0,0 +1,11 @@ +#include "drift_diffusion_updater.hpp" + +namespace bart::formulation::updater { + +template +DriftDiffusionUpdater::DriftDiffusionUpdater(std::unique_ptr diffusion_formulation_ptr, + std::shared_ptr stamper_ptr, + std::unordered_set reflective_boundaries) + : DiffusionUpdater(std::move(diffusion_formulation_ptr), stamper_ptr, reflective_boundaries) {} + +} // 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..0f7be35ce --- /dev/null +++ b/src/formulation/updater/drift_diffusion_updater.hpp @@ -0,0 +1,25 @@ +#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" + +namespace bart::formulation::updater { + +template +class DriftDiffusionUpdater : public DiffusionUpdater { + public: + using DiffusionFormulation = typename DiffusionUpdater::DiffusionFormulationType; + using Stamper = typename DiffusionUpdater::StamperType; + + DriftDiffusionUpdater(std::unique_ptr, + std::shared_ptr, + std::unordered_set reflective_boundaries = {}); + virtual ~DriftDiffusionUpdater() = default; +}; + +} // namespace bart::formulation::updater + +#endif //BART_SRC_FORMULATION_UPDATER_DRIFT_DIFFUSION_UPDATER_HPP_ 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..ecb236f95 --- /dev/null +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -0,0 +1,25 @@ +#include "formulation/updater/drift_diffusion_updater.hpp" + +#include "formulation/updater/tests/updater_tests.h" +#include "test_helpers/gmock_wrapper.h" + +namespace { + +using namespace bart; +template +using UpdaterTest = bart::formulation::updater::test_helpers::UpdaterTests; + + +template +class FormulationUpdaterDriftDiffusionTest : public UpdaterTest { + public: + static constexpr int dim{ DimensionWrapper::value }; +}; + +TYPED_TEST_SUITE(FormulationUpdaterDriftDiffusionTest, bart::testing::AllDimensions); + +TYPED_TEST(FormulationUpdaterDriftDiffusionTest, Dummy) { + EXPECT_TRUE(false); +} + +} // namespace From 32a46c1c2f441563320fe6365d07e9fe84f7aebf Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 14 Dec 2020 09:54:43 -0800 Subject: [PATCH 18/23] Added mock drift diffusion updater. --- .../scalar/tests/drift_diffusion_mock.hpp | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 src/formulation/scalar/tests/drift_diffusion_mock.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_ From 4b35d091364ef3e92e7d5a1eb5604a8c54a57cb5 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 14 Dec 2020 16:14:29 -0800 Subject: [PATCH 19/23] Added dependencies to DriftDiffusionUpdater. --- .../updater/drift_diffusion_updater.cpp | 23 +++- .../updater/drift_diffusion_updater.hpp | 21 +++- .../tests/drift_diffusion_updater_test.cpp | 100 +++++++++++++++++- 3 files changed, 137 insertions(+), 7 deletions(-) diff --git a/src/formulation/updater/drift_diffusion_updater.cpp b/src/formulation/updater/drift_diffusion_updater.cpp index 7ce72789f..4726ef211 100644 --- a/src/formulation/updater/drift_diffusion_updater.cpp +++ b/src/formulation/updater/drift_diffusion_updater.cpp @@ -3,9 +3,24 @@ namespace bart::formulation::updater { template -DriftDiffusionUpdater::DriftDiffusionUpdater(std::unique_ptr diffusion_formulation_ptr, - std::shared_ptr stamper_ptr, - std::unordered_set reflective_boundaries) - : DiffusionUpdater(std::move(diffusion_formulation_ptr), stamper_ptr, reflective_boundaries) {} +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, + 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), + 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); +} + +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 index 0f7be35ce..703dd68f0 100644 --- a/src/formulation/updater/drift_diffusion_updater.hpp +++ b/src/formulation/updater/drift_diffusion_updater.hpp @@ -5,19 +5,38 @@ #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 "utility/has_dependencies.h" namespace bart::formulation::updater { template -class DriftDiffusionUpdater : public DiffusionUpdater { +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 Stamper = typename DiffusionUpdater::StamperType; DriftDiffusionUpdater(std::unique_ptr, + std::unique_ptr, std::shared_ptr, + std::unique_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 integrated_flux_calculator_ptr() const -> IntegratedFluxCalculator* { + return integrated_flux_calculator_ptr_.get(); } + protected: + AngularFluxStorageMap angular_flux_storage_map_{}; + std::unique_ptr drift_diffusion_formulation_ptr_{ nullptr }; + std::unique_ptr integrated_flux_calculator_ptr_{ nullptr }; }; } // namespace bart::formulation::updater diff --git a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp index ecb236f95..13724e21f 100644 --- a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -1,7 +1,15 @@ #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 "system/solution/solution_types.h" namespace { @@ -10,16 +18,104 @@ template using UpdaterTest = bart::formulation::updater::test_helpers::UpdaterTests; +using ::testing::ContainerEq; + 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; + + // Supporting objects + AngularFluxStorageMap angular_flux_storage_map_{}; + + // 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 }; + + // Test parameters (each should be a different value to ensure the correct values are being passed + const int n_groups_{ test_helpers::RandomInt(1, 4) }; + const Vector::size_type angular_flux_size{ static_cast(test_helpers::RandomInt(5, 10)) }; + const int n_angles_{ n_groups_ + 1 }; + + 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::make_unique(); + stamper_obs_ptr_ = stamper_ptr.get(); + + using Index = system::SolutionIndex; + for (int group = 0; group < n_groups_; ++group) { + for (int angle = 0; angle < n_angles_; ++angle) { + Index solution_index{ system::EnergyGroup(group), system::AngleIdx(angle) }; + angular_flux_storage_map_.insert({solution_index, std::make_shared(angular_flux_size)}); + } + } +} + TYPED_TEST_SUITE(FormulationUpdaterDriftDiffusionTest, bart::testing::AllDimensions); -TYPED_TEST(FormulationUpdaterDriftDiffusionTest, Dummy) { - EXPECT_TRUE(false); +// 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_unique(), + std::make_unique(), + 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_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; + + std::unique_ptr test_updater{ nullptr }; + const int n_dependencies_to_check{ 4 }; + 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(), + this->angular_flux_storage_map_); + }); + } } } // namespace From fb9d1b38c43482b07d4712da347e1e7ab80fccdc Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 15 Dec 2020 10:11:54 -0800 Subject: [PATCH 20/23] Added initial implementation of SetUpFixedFunctions to DriftDiffusionUpdater. --- src/formulation/updater/diffusion_updater.cpp | 1 - .../updater/drift_diffusion_updater.cpp | 6 ++ .../updater/drift_diffusion_updater.hpp | 2 + .../tests/drift_diffusion_updater_test.cpp | 57 ++++++++++++++++++- 4 files changed, 62 insertions(+), 4 deletions(-) diff --git a/src/formulation/updater/diffusion_updater.cpp b/src/formulation/updater/diffusion_updater.cpp index 6ccfd91e2..b2bef1ec9 100644 --- a/src/formulation/updater/diffusion_updater.cpp +++ b/src/formulation/updater/diffusion_updater.cpp @@ -34,7 +34,6 @@ template 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 { diff --git a/src/formulation/updater/drift_diffusion_updater.cpp b/src/formulation/updater/drift_diffusion_updater.cpp index 4726ef211..5aadb550d 100644 --- a/src/formulation/updater/drift_diffusion_updater.cpp +++ b/src/formulation/updater/drift_diffusion_updater.cpp @@ -18,6 +18,12 @@ DriftDiffusionUpdater::DriftDiffusionUpdater( AssertPointerNotNull(drift_diffusion_formulation_ptr_.get(), "drift diffusion formulation", call_location); AssertPointerNotNull(integrated_flux_calculator_ptr_.get(), "integrated flux calculator", 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); +} template class DriftDiffusionUpdater<1>; template class DriftDiffusionUpdater<2>; diff --git a/src/formulation/updater/drift_diffusion_updater.hpp b/src/formulation/updater/drift_diffusion_updater.hpp index 703dd68f0..a80f5115d 100644 --- a/src/formulation/updater/drift_diffusion_updater.hpp +++ b/src/formulation/updater/drift_diffusion_updater.hpp @@ -27,6 +27,7 @@ class DriftDiffusionUpdater : public DiffusionUpdater, public utility::HasD 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* { @@ -34,6 +35,7 @@ class DriftDiffusionUpdater : public DiffusionUpdater, public utility::HasD 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::unique_ptr drift_diffusion_formulation_ptr_{ nullptr }; std::unique_ptr integrated_flux_calculator_ptr_{ nullptr }; diff --git a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp index 13724e21f..98d6a98fe 100644 --- a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -9,6 +9,7 @@ #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" namespace { @@ -18,7 +19,7 @@ template using UpdaterTest = bart::formulation::updater::test_helpers::UpdaterTests; -using ::testing::ContainerEq; +using ::testing::ContainerEq, ::testing::DoDefault, ::testing::_, ::testing::Ref; template class FormulationUpdaterDriftDiffusionTest : public UpdaterTest { @@ -32,6 +33,10 @@ class FormulationUpdaterDriftDiffusionTest : public UpdaterTest; + // Test object + using TestUpdater = formulation::updater::DriftDiffusionUpdater; + std::unique_ptr test_updater_ptr_{ nullptr }; + // Supporting objects AngularFluxStorageMap angular_flux_storage_map_{}; @@ -60,7 +65,7 @@ void FormulationUpdaterDriftDiffusionTest::SetUp() { 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::make_unique(); + auto stamper_ptr = std::shared_ptr(this->MakeStamper()); stamper_obs_ptr_ = stamper_ptr.get(); using Index = system::SolutionIndex; @@ -70,6 +75,12 @@ void FormulationUpdaterDriftDiffusionTest::SetUp() { angular_flux_storage_map_.insert({solution_index, std::make_shared(angular_flux_size)}); } } + 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), + angular_flux_storage_map_, + reflective_boundaries); } TYPED_TEST_SUITE(FormulationUpdaterDriftDiffusionTest, bart::testing::AllDimensions); @@ -87,7 +98,7 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorDependencyGetters) { EXPECT_NO_THROW({ test_updater = std::make_unique(std::make_unique(), std::make_unique(), - std::make_unique(), + std::make_shared(), std::make_unique(), this->angular_flux_storage_map_); }); @@ -118,4 +129,44 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorBadDependencies) { } } +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()); + + 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)); + 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(2) + .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 From cd56dbeadaf04115d005df4f61e428108cbd6fc9 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Tue, 15 Dec 2020 11:01:10 -0800 Subject: [PATCH 21/23] Completed implementation of DriftDiffusionUpdater. --- .../updater/drift_diffusion_updater.cpp | 19 ++++++++ .../updater/drift_diffusion_updater.hpp | 6 +++ .../tests/drift_diffusion_updater_test.cpp | 47 +++++++++++++++---- 3 files changed, 63 insertions(+), 9 deletions(-) diff --git a/src/formulation/updater/drift_diffusion_updater.cpp b/src/formulation/updater/drift_diffusion_updater.cpp index 5aadb550d..dbbcb63dd 100644 --- a/src/formulation/updater/drift_diffusion_updater.cpp +++ b/src/formulation/updater/drift_diffusion_updater.cpp @@ -8,21 +8,40 @@ DriftDiffusionUpdater::DriftDiffusionUpdater( 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>; diff --git a/src/formulation/updater/drift_diffusion_updater.hpp b/src/formulation/updater/drift_diffusion_updater.hpp index a80f5115d..846fb5b8f 100644 --- a/src/formulation/updater/drift_diffusion_updater.hpp +++ b/src/formulation/updater/drift_diffusion_updater.hpp @@ -6,6 +6,7 @@ #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 { @@ -17,12 +18,14 @@ class DriftDiffusionUpdater : public DiffusionUpdater, public utility::HasD 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; @@ -32,11 +35,14 @@ class DriftDiffusionUpdater : public DiffusionUpdater, public utility::HasD 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 }; }; diff --git a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp index 98d6a98fe..ed2743c92 100644 --- a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -11,6 +11,7 @@ #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 { @@ -19,7 +20,7 @@ template using UpdaterTest = bart::formulation::updater::test_helpers::UpdaterTests; -using ::testing::ContainerEq, ::testing::DoDefault, ::testing::_, ::testing::Ref; +using ::testing::ContainerEq, ::testing::DoDefault, ::testing::_, ::testing::Ref, ::testing::Return, ::testing::ReturnRef; template class FormulationUpdaterDriftDiffusionTest : public UpdaterTest { @@ -32,24 +33,28 @@ class FormulationUpdaterDriftDiffusionTest : public UpdaterTest; 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_{}; + 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 int n_groups_{ test_helpers::RandomInt(1, 4) }; const Vector::size_type angular_flux_size{ static_cast(test_helpers::RandomInt(5, 10)) }; - const int n_angles_{ n_groups_ + 1 }; std::unordered_set reflective_boundaries{Boundary::kXMin, Boundary::kYMax}; @@ -67,18 +72,27 @@ void FormulationUpdaterDriftDiffusionTest::SetUp() { 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 < n_groups_; ++group) { - for (int angle = 0; angle < n_angles_; ++angle) { + 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, std::make_shared(angular_flux_size)}); + 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); } @@ -100,11 +114,13 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorDependencyGetters) { 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); + 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_)); } @@ -115,15 +131,18 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorBadDependencies) { 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{ 4 }; + 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_); }); } @@ -136,11 +155,21 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, UpdateFixedTermTest) { 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) { @@ -159,7 +188,7 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, UpdateFixedTermTest) { } EXPECT_CALL(*this->stamper_obs_ptr_, StampMatrix(Ref(*this->matrix_to_stamp), _)) - .Times(2) + .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()); From 7066f5c20bf18b617db6f062c8bc7a60b08afd84 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 4 Jan 2021 22:42:48 -0800 Subject: [PATCH 22/23] added missing test to verify drift_diffusion_updater dependency --- src/formulation/updater/tests/drift_diffusion_updater_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp index ed2743c92..74b7c5a0b 100644 --- a/src/formulation/updater/tests/drift_diffusion_updater_test.cpp +++ b/src/formulation/updater/tests/drift_diffusion_updater_test.cpp @@ -120,6 +120,7 @@ TYPED_TEST(FormulationUpdaterDriftDiffusionTest, ConstructorDependencyGetters) { 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_)); } From 070fa9e8bd33e4dace64c8cf17e29d2eee0d6d53 Mon Sep 17 00:00:00 2001 From: Joshua Rehak Date: Mon, 4 Jan 2021 22:47:36 -0800 Subject: [PATCH 23/23] added coverage exclusion for part of fixed updater not used by any formulations --- src/formulation/updater/fixed_updater.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/formulation/updater/fixed_updater.hpp b/src/formulation/updater/fixed_updater.hpp index aead90bf3..faec14751 100644 --- a/src/formulation/updater/fixed_updater.hpp +++ b/src/formulation/updater/fixed_updater.hpp @@ -36,8 +36,11 @@ class FixedUpdater : public FixedUpdaterI { 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: