Feature update v0.1.0, SAAF Refl BC, 1D Ang. Quad, Fixed Source Solve

.gitignore

@@ -34,3 +34,4 @@ html/
 latex/
 build/
 test_data/
+*.csv

.travis.yml

@@ -5,26 +5,30 @@ notifications:
 
 language: cpp
 
+env:
+ global:
+  - DOCKER_TAG=v9.1.1
+
 services:
  - docker
 
 before_install:
- - sudo docker pull jsrehak/bart
+ - sudo docker pull jsrehak/bart:$DOCKER_TAG
  - ci_env=`bash <(curl -s https://codecov.io/env)`
 
 install:
- - sudo docker run jsrehak/bart
+ - sudo docker run jsrehak/bart:$DOCKER_TAG
  - id=$(sudo docker ps -l | tail -n1 | awk '{print $1}')
  - sudo docker cp ../BART $id:/home/bart
  - id=$(sudo docker ps -l | tail -n1 | awk '{print $1}')
- - sudo docker commit $id jsrehak/bart
- - sudo docker run $ci_env -it -u="root" -d -w="/home/bart/BART" --name bart jsrehak/bart
+ - sudo docker commit $id jsrehak/bart:$DOCKER_TAG
+ - sudo docker run $ci_env -it -u="root" -d -w="/home/bart/BART" --name bart jsrehak/bart:$DOCKER_TAG
 
 script:
  - docker exec bart cmake .
  - docker exec bart make -j4
  - docker exec bart bash -c "./bart_test"
- - docker exec bart bash -c "mpirun -np 2 --allow-run-as-root ./bart_test --mpi -l 0"
+ - docker exec bart bash -c "mpirun -np 2 --allow-run-as-root --oversubscribe ./bart_test --mpi -l 0"
 
 after_success:
  - docker exec bart bash -c "./coverage.sh"

README.md

@@ -36,8 +36,8 @@ BART uses MPI for parallelism and is designed for computation on distributed mem
 ### Formulations & Methods
 
 BART supports the following formulations: 
-- the diffusion equation in 1/2/3D.
-- the Self-Adjoint Angular Flux formulation in 1/2/3D (Only 3D is currently usable because only a 3D angular quadrature set has been implemented).
+- the diffusion equation in 1/2/3D
+- the Self-Adjoint Angular Flux formulation in 1D & 3D (2D pending implementation of a 2D quadrature)
 
 One of the major design goals of BART is to provide a framework for testing acceleration methods. There are no methods implemented in the current version, but multiple methods are planned to be implemented, including:
 

benchmarks/picca_2016/figure_2.material_map

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

benchmarks/picca_2016/figure_2_diffusion.input

@@ -0,0 +1,20 @@
+set problem dimension                        = 1
+set transport model                          = diffusion
+set number of groups                         = 1
+set angular quadrature name                  = gauss_legendre
+set angular quadrature order                 = 2
+set do eigenvalue calculations               = true
+
+set x, y, z max values of boundary locations = 100.0
+set number of cells for x, y, z directions   = 20
+set uniform refinements                      = 2
+set number of materials                      = 2
+
+set finite element polynomial degree         = 1
+
+set output file name base                    = figure_2_diffusion
+
+subsection material ID map
+set material id file name map                = 1: reactor, 2: reflector
+set material id file name                    = figure_2.material_map
+end

benchmarks/picca_2016/figure_2_diffusion_half.input

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

benchmarks/picca_2016/figure_2_half.material_map

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

benchmarks/picca_2016/figure_2_saaf.input

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

benchmarks/picca_2016/figure_2_saaf_half.input

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

benchmarks/sood_1999/pu_2_0_in/pu_2_0_in.input

@@ -1,13 +1,13 @@
-set problem dimension                        = 2
+set problem dimension                        = 1
 set transport model                          = diffusion
 set number of groups                         = 2
 set do eigenvalue calculations               = true
-set have reflective boundary                 = true
-set reflective boundary names                = xmin, xmax, ymax, ymin
+set reflective boundary names                = xmin, xmax
 
-set x, y, z max values of boundary locations = 10.0, 10.0
-set number of cells for x, y, z directions   = 10, 10
+set x, y, z max values of boundary locations = 10.0
+set number of cells for x, y, z directions   = 10
 set number of materials                      = 1
+set uniform refinements                      = 2
 
 set finite element polynomial degree         = 1
 

benchmarks/sood_1999/pu_2_0_in_saaf/pu_2_0_in.input → benchmarks/sood_1999/pu_2_0_in_saaf/pu_2_0_in_saaf.input

@@ -1,19 +1,20 @@
-set problem dimension                        = 3
+set problem dimension                        = 1
 set transport model                          = saaf
 set number of groups                         = 2
-set angular quadrature name                  = level_symmetric_gaussian
+set angular quadrature name                  = gauss_legendre
 set angular quadrature order                 = 2
 set do eigenvalue calculations               = true
 set have reflective boundary                 = true
-set reflective boundary names                = xmin, xmax, ymax, ymin
+set reflective boundary names                = xmin, xmax
 
-set x, y, z max values of boundary locations = 10.0, 10.0, 10.0
-set number of cells for x, y, z directions   = 10, 10, 10
+set x, y, z max values of boundary locations = 100.0
+set number of cells for x, y, z directions   = 10
+set uniform refinements                      = 1
 set number of materials                      = 1
 
 set finite element polynomial degree         = 1
 
-set output file name base                    = pu_2_0_in
+set output file name base                    = pu_2_0_in_saaf
 
 subsection material ID map
 set material id file name map                = 1: pu_2_0_in_pu_239.material

src/convergence/final.cc

@@ -30,6 +30,7 @@ Final<CompareType>& Final<CompareType>::SetIteration(
 
 template class Final<system::moments::MomentVector>;
 template class Final<system::moments::MomentsMap>;
+template class Final<const system::moments::MomentsMap>;
 template class Final<double>;
 
 

src/convergence/final_checker_or_n.cc

@@ -55,6 +55,8 @@ template class FinalCheckerOrN<system::moments::MomentVector,
                                moments::SingleMomentCheckerI>;
 template class FinalCheckerOrN<system::moments::MomentsMap,
                                moments::MultiMomentCheckerI>;
+template class FinalCheckerOrN<const system::moments::MomentsMap,
+                               moments::MultiMomentCheckerI>;
 template class FinalCheckerOrN<double, parameters::SingleParameterChecker>;
 
 

src/convergence/moments/single_moment_checker_l1_norm.cc

@@ -9,7 +9,7 @@ namespace moments {
 bool SingleMomentCheckerL1Norm::CheckIfConverged(
     const system::moments::MomentVector &current_iteration,
     const system::moments::MomentVector &previous_iteration) {
-  system::moments::MomentVector difference{current_iteration};
+  system::moments::MomentVector difference(current_iteration);
   difference -= previous_iteration;
   delta_ = difference.l1_norm()/current_iteration.l1_norm();
   is_converged_ = delta_ <= max_delta_;

src/domain/definition.cc

@@ -55,11 +55,17 @@ Definition<1>::Definition(
 template <int dim>
 Definition<dim>& Definition<dim>::SetUpMesh() {
   // SetUp triangulation using the mesh dependency
+  return SetUpMesh(0);
+}
+
+template<int dim>
+Definition<dim>& Definition<dim>::SetUpMesh(const int global_refinements) {
   AssertThrow(mesh_->has_material_mapping(),
-                    dealii::ExcMessage("Mesh object must have initialized material mapping"));
+              dealii::ExcMessage("Mesh object must have initialized material mapping"));
   mesh_->FillTriangulation(triangulation_);
   mesh_->FillBoundaryID(triangulation_);
   mesh_->FillMaterialID(triangulation_);
+  triangulation_.refine_global(global_refinements);
   return *this;
 }
 
@@ -170,7 +176,6 @@ int Definition<dim>::total_degrees_of_freedom() const {
   return total_degrees_of_freedom_;
 }
 
-
 template class Definition<1>;
 template class Definition<2>;
 template class Definition<3>;

src/domain/definition.h

@@ -7,7 +7,7 @@
 #include <deal.II/base/index_set.h>
 #include <deal.II/base/iterator_range.h>
 #include <deal.II/distributed/tria.h>
-#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
 
 #include "domain/definition_i.h"
@@ -74,6 +74,7 @@ class Definition : public DefinitionI<dim> {
 
   Definition<dim>& SetUpDOF() override;
   Definition<dim>& SetUpMesh() override;
+  Definition<dim>& SetUpMesh(const int global_refinements) override;
 
   dealii::FullMatrix<double> GetCellMatrix() const override {
     int cell_dofs = finite_element_->dofs_per_cell();
@@ -132,7 +133,7 @@ class Definition : public DefinitionI<dim> {
   dealii::IndexSet locally_relevant_dofs_;
 
   /*! Constraint matrix */
-  dealii::ConstraintMatrix constraint_matrix_;
+  dealii::AffineConstraints<double> constraint_matrix_;
 
   /*! Dynamic sparsity pattern for MPI matrices */
   dealii::DynamicSparsityPattern dynamic_sparsity_pattern_;

src/domain/definition_i.h

@@ -47,12 +47,21 @@ class DefinitionI : public utility::HasDescription {
   /*! Set up the DOF handler, to access sparsity patterns, etc */
   virtual DefinitionI<dim>& SetUpDOF() = 0;
 
-/*! Fills triangulation with mesh defined in MeshI object
+  /*! \brief Fills triangulation with mesh defined in MeshI object
    * Creates mesh shape, sets up boundary ids and material ids. Requires that
    * the mesh has a material mapping setup.
    */
   virtual DefinitionI<dim>& SetUpMesh() = 0;
 
+  /*! \brief Fills triangulation with mesh defined in MeshI object.
+   *
+   * Creates mesh shape, sets up boundary ids and material ids. Requires that
+   * the mesh has a material mapping setup. Also performs global refinements.
+   *
+   * @param global_refinements number of global refinements to perform
+   */
+  virtual DefinitionI<dim>& SetUpMesh(const int global_refinements) = 0;
+
   /*! Get a matrix suitible for a cell matrix.
    *
    * \return a dealii FullMatrix<double> of appropriate size.

src/domain/finite_element/finite_element.cc

@@ -47,13 +47,21 @@ template<int dim>
 std::vector<double> FiniteElement<dim>::ValueAtQuadrature(
     const system::moments::MomentVector moment) const {
 
-  std::vector<double> return_vector(finite_element_->dofs_per_cell, 0);
+  std::vector<double> return_vector(n_cell_quad_pts(), 0);
 
   values_->get_function_values(moment, return_vector);
 
   return return_vector;
 }
 
+template<int dim>
+std::vector<double> FiniteElement<dim>::ValueAtFaceQuadrature(
+    const dealii::Vector<double>& values_at_dofs) const {
+  std::vector<double> return_vector(n_face_quad_pts(), 0);
+  face_values_->get_function_values(values_at_dofs, return_vector);
+  return return_vector;
+}
+
 template class FiniteElement<1>;
 template class FiniteElement<2>;
 template class FiniteElement<3>;

src/domain/finite_element/finite_element.h

@@ -3,6 +3,7 @@
 
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/fe/fe_values.h>
+#include <system/system_types.h>
 
 #include "domain/finite_element/finite_element_i.h"
 
@@ -72,7 +73,11 @@ class FiniteElement : public FiniteElementI<dim> {
     return face_values_->normal_vector(0);
   };
 
-  std::vector<double> ValueAtQuadrature(const system::moments::MomentVector moment) const override;
+  std::vector<double> ValueAtQuadrature(
+      const system::moments::MomentVector moment) const override;
+
+  std::vector<double> ValueAtFaceQuadrature(
+      const dealii::Vector<double>& values_at_dofs) const override;
 
  protected:
   std::shared_ptr<dealii::FiniteElement<dim, dim>> finite_element_;

src/domain/finite_element/finite_element_i.h

@@ -5,6 +5,7 @@
 
 #include "domain/domain_types.h"
 #include "system/moments/spherical_harmonic_types.h"
+#include "system/system_types.h"
 
 /*! \brief Interface for a finite element object based on the dealii library.
  *
@@ -137,6 +138,14 @@ class FiniteElementI {
   virtual std::vector<double> ValueAtQuadrature(
       const system::moments::MomentVector moment) const = 0;
 
+  /*! \brief Get the value of an MPI Vector at the cell face quadrature points.
+   *
+   * @param mpi_vector mpi vector to get the face values of.
+   * @return a vector holding the value of the mpi vector at each face quadrature point.
+   */
+   virtual std::vector<double> ValueAtFaceQuadrature(
+       const dealii::Vector<double>& values_at_dofs) const = 0;
+
   // DealII Finite element object access. These methods access the underlying
   // finite element objects.
   /*! \brief Gets pointer to the underlying finite element object */

src/domain/finite_element/tests/finite_element_gaussian_test.cc

@@ -200,4 +200,10 @@ TYPED_TEST(DomainFiniteElementGaussianBaseMethodsTest, BaseValueAtQuadrature) {
   this->TestValueAtQuadrature(&test_fe);
 }
 
+TYPED_TEST(DomainFiniteElementGaussianBaseMethodsTest, BaseValueAtFaceQuadrature) {
+  bart::domain::finite_element::FiniteElementGaussian<this->dim> test_fe{
+      problem::DiscretizationType::kDiscontinuousFEM, 2};
+  this->TestValueAtFaceQuadrature(&test_fe);
+}
+
 } // namespace

src/domain/finite_element/tests/finite_element_mock.h

@@ -45,6 +45,9 @@ class FiniteElementMock : public FiniteElementI<dim> {
 
   MOCK_METHOD(std::vector<double>, ValueAtQuadrature, (const system::moments::MomentVector moment), (const, override));
 
+  MOCK_METHOD(std::vector<double>, ValueAtFaceQuadrature,
+      (const dealii::Vector<double>&), (const, override));
+
   MOCK_METHOD((dealii::FiniteElement<dim, dim>*), finite_element, (), (override));
 
   MOCK_METHOD(dealii::FEValues<dim>*, values, (), (override));