Merge pull request #1 from satfra/discretization_renaming

Discretization renaming

CHANGELOG.md

@@ -0,0 +1,12 @@
+# Changelog
+
+## Version 1.0.1
+
+### Fixed
+
+- There was a bug with TaskFlow internally in deal.ii. Fixed for now by simply disabling taskflow in deal.ii.
+- deal.ii changed its interface for dealii::SolutionTransfer. Adapted the corresponding methods.
+
+### Changed
+
+- The FlowingVariables classes are now in separate namespaces. For finite elements, use DiFfRG::FE::FlowingVariables, for pure variable systems use DiFfRG::FlowingVariables.

DiFfRG/documentation/tutorials/tut1.md

@@ -137,7 +137,7 @@ We now use the types we defined above to construct objects of all the classes de
 With everything prepared, we are now ready to set up and run the equation system. To do so, we create an initial condition on the discretization (i.e. FEM space) we set up before,
 ```Cpp
   // Set up the initial condition
-  FlowingVariables initial_condition(discretization);
+  FE::FlowingVariables initial_condition(discretization);
   initial_condition.interpolate(model);
 ```
 and use it to run the time-stepper from RG-time 0 to the final RG-time, which we infer from the parameter file:

DiFfRG/include/DiFfRG/discretization/FEM/cg.hh

@@ -80,6 +80,13 @@ namespace DiFfRG
 
       void reinit() { setup_dofs(); }
 
+      std::vector<uint> get_block_structure() const
+      {
+        std::vector<uint> block_structure{dof_handler.n_dofs()};
+        if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables());
+        return block_structure;
+      }
+
     protected:
       void setup_dofs()
       {

DiFfRG/include/DiFfRG/discretization/FEM/dg.hh

@@ -94,6 +94,13 @@ namespace DiFfRG
         return dof;
       }
 
+      std::vector<uint> get_block_structure() const
+      {
+        std::vector<uint> block_structure{dof_handler.n_dofs()};
+        if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables());
+        return block_structure;
+      }
+
     protected:
       void setup_dofs()
       {

DiFfRG/include/DiFfRG/discretization/FEM/ldg.hh

@@ -88,6 +88,13 @@ namespace DiFfRG
         return dof;
       }
 
+      std::vector<uint> get_block_structure() const
+      {
+        std::vector<uint> block_structure{dof_handler[0]->n_dofs()};
+        if (Components::count_variables() > 0) block_structure.push_back(Components::count_variables());
+        return block_structure;
+      }
+
     protected:
       void setup_dofs()
       {

DiFfRG/include/DiFfRG/discretization/FV/discretization.hh

@@ -0,0 +1,56 @@
+#pragma once
+
+// external libraries
+#include <deal.II/base/point.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/lac/affine_constraints.h>
+#include <spdlog/spdlog.h>
+
+// DiFfRG
+#include <DiFfRG/common/utils.hh>
+
+namespace DiFfRG
+{
+  namespace FV
+  {
+    using namespace dealii;
+
+    /**
+     * @brief Class to manage the system on which we solve, i.e. fe spaces, grids, etc.
+     * This class is a System for CG systems.
+     *
+     * @tparam Model_ The Model class used for the Simulation
+     */
+    template <typename Components_, typename NumberType_, typename Mesh_> class Discretization
+    {
+    public:
+      using Components = Components_;
+      using NumberType = NumberType_;
+      using VectorType = Vector<NumberType>;
+      using SparseMatrixType = SparseMatrix<NumberType>;
+      using Mesh = Mesh_;
+      static constexpr uint dim = Mesh::dim;
+
+      Discretization(Mesh &mesh, const JSONValue &json) : mesh(mesh), json(json) {};
+
+      const auto &get_mapping() const { return mapping; }
+      const auto &get_triangulation() const { return mesh.get_triangulation(); }
+      auto &get_triangulation() { return mesh.get_triangulation(); }
+      const auto &get_json() const { return json; }
+
+      void reinit() {}
+
+    protected:
+      Mesh &mesh;
+      JSONValue json;
+
+      AffineConstraints<NumberType> constraints;
+      MappingQ1<dim> mapping;
+    };
+  } // namespace FV
+} // namespace DiFfRG

DiFfRG/include/DiFfRG/discretization/data/data.hh

@@ -27,21 +27,6 @@ namespace DiFfRG
     private:
       FUN fun;
     };
-
-    struct NoDiscretization {
-      static constexpr uint dim = 0;
-      using NumberType = double;
-      using Components = void;
-
-      const DoFHandler<dim> &get_dof_handler() const { throw std::runtime_error("No discretization chosen!"); }
-
-      template <typename Model> static std::vector<uint> get_block_structure(const Model &)
-      {
-        std::vector<uint> block_structure{0};
-        block_structure.push_back(Model::Components::count_variables());
-        return block_structure;
-      }
-    };
   } // namespace internal
 
   /**
@@ -51,29 +36,15 @@ namespace DiFfRG
    *
    * @tparam Discretization Spatial Discretization used in the system
    */
-  template <typename Discretization = typename internal::NoDiscretization>
-  class FlowingVariables : public AbstractFlowingVariables<typename Discretization::NumberType>
+  template <typename NT = double> class FlowingVariables : public AbstractFlowingVariables<NT>
   {
   public:
-    using NumberType = typename Discretization::NumberType;
-    using Components = typename Discretization::Components;
-    static constexpr uint dim = Discretization::dim;
+    using NumberType = NT;
 
     /**
      * @brief Construct a new Flowing Variables object
-     *
-     * @param discretization The spatial discretization to use
      */
-    FlowingVariables(const Discretization &discretization) : dof_handler(discretization.get_dof_handler())
-    {
-      auto block_structure = Components::get_block_structure(dof_handler);
-      m_data = (block_structure);
-    }
-
-    /**
-     * @brief Construct a new Flowing Variables object
-     */
-    FlowingVariables() : dof_handler(*static_cast<DoFHandler<dim> *>(0)) {}
+    FlowingVariables() {}
 
     /**
      * @brief Interpolates the initial condition from a numerical model.
@@ -83,15 +54,9 @@ namespace DiFfRG
      */
     template <typename Model> void interpolate(const Model &model)
     {
-      if constexpr (Model::Components::count_fe_functions() > 0) {
-        auto interpolating_function = [&model](const auto &p, auto &values) { model.initial_condition(p, values); };
-        internal::FunctionFromLambda<dim, NumberType> initial_condition_function(std::move(interpolating_function),
-                                                                                 dof_handler.get_fe().n_components());
-        VectorTools::interpolate(dof_handler, initial_condition_function, m_data.block(0));
-      } else {
-        auto block_structure = Discretization::get_block_structure(model);
-        m_data = (block_structure);
-      }
+      std::vector<uint> block_structure{0};
+      block_structure.push_back(Model::Components::count_variables());
+      m_data = (block_structure);
       if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1));
     }
 
@@ -120,8 +85,156 @@ namespace DiFfRG
     virtual const Vector<NumberType> &variable_data() const override { return m_data.block(1); }
 
   private:
-    const DoFHandler<dim> &dof_handler;
     BlockVector<NumberType> m_data;
   };
 
+  namespace FE
+  {
+    /**
+     * @brief A class to set up initial data for whatever discretization we have chosen.
+     *        Also used to switch/manage memory, vectors, matrices over interfaces between spatial discretization and
+     * separate variables.
+     *
+     * @tparam Discretization Spatial Discretization used in the system
+     */
+    template <typename Discretization>
+    class FlowingVariables : public AbstractFlowingVariables<typename Discretization::NumberType>
+    {
+    public:
+      using NumberType = typename Discretization::NumberType;
+      using Components = typename Discretization::Components;
+      static constexpr uint dim = Discretization::dim;
+
+      /**
+       * @brief Construct a new Flowing Variables object
+       *
+       * @param discretization The spatial discretization to use
+       */
+      FlowingVariables(const Discretization &discretization)
+          : discretization(discretization), dof_handler(discretization.get_dof_handler())
+      {
+      }
+
+      /**
+       * @brief Interpolates the initial condition from a numerical model.
+       *
+       * @param model The model to interpolate from. Must provide a method initial_condition(const Point<dim> &,
+       * Vector<NumberType> &)
+       */
+      template <typename Model> void interpolate(const Model &model)
+      {
+        auto block_structure = discretization.get_block_structure();
+        m_data = (block_structure);
+
+        if constexpr (Model::Components::count_fe_functions() > 0) {
+          auto interpolating_function = [&model](const auto &p, auto &values) { model.initial_condition(p, values); };
+          internal::FunctionFromLambda<dim, NumberType> initial_condition_function(std::move(interpolating_function),
+                                                                                   dof_handler.get_fe().n_components());
+          VectorTools::interpolate(dof_handler, initial_condition_function, m_data.block(0));
+        }
+        if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1));
+      }
+
+      /**
+       * @brief Obtain the data vector holding both spatial (block 0) and variable (block 1) data.
+       *
+       * @return BlockVector<NumberType>& The data vector.
+       */
+      virtual BlockVector<NumberType> &data() override { return m_data; }
+      virtual const BlockVector<NumberType> &data() const override { return m_data; }
+
+      /**
+       * @brief Obtain the spatial data vector.
+       *
+       * @return Vector<NumberType>& The spatial data vector.
+       */
+      virtual Vector<NumberType> &spatial_data() override { return m_data.block(0); }
+      virtual const Vector<NumberType> &spatial_data() const override { return m_data.block(0); }
+
+      /**
+       * @brief Obtain the variable data vector.
+       *
+       * @return Vector<NumberType>& The variable data vector.
+       */
+      virtual Vector<NumberType> &variable_data() override { return m_data.block(1); }
+      virtual const Vector<NumberType> &variable_data() const override { return m_data.block(1); }
+
+    private:
+      const Discretization &discretization;
+      const DoFHandler<dim> &dof_handler;
+      BlockVector<NumberType> m_data;
+    };
+  } // namespace FE
+
+  namespace FV
+  {
+    /**
+     * @brief A class to set up initial data for whatever discretization we have chosen.
+     *        Also used to switch/manage memory, vectors, matrices over interfaces between spatial discretization and
+     * separate variables.
+     *
+     * @tparam Discretization Spatial Discretization used in the system
+     */
+    template <typename Discretization>
+    class FlowingVariables : public AbstractFlowingVariables<typename Discretization::NumberType>
+    {
+    public:
+      using NumberType = typename Discretization::NumberType;
+      using Components = typename Discretization::Components;
+      static constexpr uint dim = Discretization::dim;
+
+      /**
+       * @brief Construct a new Flowing Variables object
+       *
+       * @param discretization The spatial discretization to use
+       */
+      FlowingVariables(const Discretization &discretization) : discretization(discretization) {}
+
+      /**
+       * @brief Interpolates the initial condition from a numerical model.
+       *
+       * @param model The model to interpolate from. Must provide a method initial_condition(const Point<dim> &,
+       * Vector<NumberType> &)
+       */
+      template <typename Model> void interpolate(const Model &model)
+      {
+        auto block_structure = discretization.get_block_structure();
+        m_data = (block_structure);
+
+        if constexpr (Model::Components::count_fe_functions() > 0) {
+          // TODO
+        }
+        if (m_data.n_blocks() > 1) model.initial_condition_variables(m_data.block(1));
+      }
+
+      /**
+       * @brief Obtain the data vector holding both spatial (block 0) and variable (block 1) data.
+       *
+       * @return BlockVector<NumberType>& The data vector.
+       */
+      virtual BlockVector<NumberType> &data() override { return m_data; }
+      virtual const BlockVector<NumberType> &data() const override { return m_data; }
+
+      /**
+       * @brief Obtain the spatial data vector.
+       *
+       * @return Vector<NumberType>& The spatial data vector.
+       */
+      virtual Vector<NumberType> &spatial_data() override { return m_data.block(0); }
+      virtual const Vector<NumberType> &spatial_data() const override { return m_data.block(0); }
+
+      /**
+       * @brief Obtain the variable data vector.
+       *
+       * @return Vector<NumberType>& The variable data vector.
+       */
+      virtual Vector<NumberType> &variable_data() override { return m_data.block(1); }
+      virtual const Vector<NumberType> &variable_data() const override { return m_data.block(1); }
+
+    private:
+      const Discretization &discretization;
+      BlockVector<NumberType> m_data;
+    };
+  } // namespace FV
+
 } // namespace DiFfRG

DiFfRG/include/DiFfRG/discretization/mesh/h_adaptivity.hh

@@ -101,7 +101,7 @@ namespace DiFfRG
       assembler.reinit();
       assembler.reinit_vector(solution);
 
-      solution_trans.interpolate(previous_solution, solution);
+      solution_trans.interpolate(solution);
       constraints.distribute(solution);
 
       return true;

DiFfRG/include/DiFfRG/model/component_descriptor.hh

@@ -165,13 +165,6 @@ namespace DiFfRG
       return dofs_per_component;
     }
 
-    template <typename DoFH> static std::vector<uint> get_block_structure(const DoFH &dofh)
-    {
-      std::vector<uint> block_structure{dofh.n_dofs()};
-      if (n_variables > 0) block_structure.push_back(n_variables);
-      return block_structure;
-    }
-
   private:
     template <typename T> using SubsystemMatrix = std::array<std::array<T, n_fe_subsystems>, n_fe_subsystems>;
 

Examples/ONfiniteT/CG.cc

@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
   TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor);
 
   // Set up the initial condition
-  FlowingVariables initial_condition(discretization);
+  FE::FlowingVariables initial_condition(discretization);
   initial_condition.interpolate(model);
 
   // Now we start the timestepping

Examples/ONfiniteT/LDG.cc

@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
   TimeStepper time_stepper(json, &assembler, &data_out, &mesh_adaptor);
 
   // Set up the initial condition
-  FlowingVariables initial_condition(discretization);
+  FE::FlowingVariables initial_condition(discretization);
   initial_condition.interpolate(model);
 
   // Now we start the timestepping