Begining of the code for unsteady heat

src/language/algorithms/CMakeLists.txt

@@ -5,6 +5,8 @@ add_library(PugsLanguageAlgorithms
   ElasticityDiamondAlgorithm.cpp
   Heat5PointsAlgorithm.cpp
   HeatDiamondAlgorithm.cpp
+  HeatDiamondAlgorithm2.cpp
+  ParabolicHeat.cpp
 )
 
 

src/language/algorithms/ParabolicHeat.cpp

+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LeastSquareSolver.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/LocalRectangularMatrix.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/algorithms/ParabolicHeat.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/ScalarDiamondScheme.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Timer.hpp>
+
+template <size_t Dimension>
+ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
+  std::shared_ptr<const IMesh> i_mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& T_id,
+  const FunctionSymbolId& kappa_id,
+  const FunctionSymbolId& f_id,
+  const double& final_time,
+  const double& dt)
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+  using MeshDataType     = MeshData<Dimension>;
+
+  constexpr ItemType FaceType = [] {
+    if constexpr (Dimension > 1) {
+      return ItemType::face;
+    } else {
+      return ItemType::node;
+    }
+  }();
+  // HERE build CLs, time loop, main outputs
+  using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition,
+                                         SymmetryBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+  BoundaryConditionList boundary_condition_list;
+
+  std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  for (const auto& bc_descriptor : bc_descriptor_list) {
+    bool is_valid_boundary_condition = true;
+
+    switch (bc_descriptor->type()) {
+    case IBoundaryConditionDescriptor::Type::symmetry: {
+      throw NotImplementedError("NIY");
+      break;
+    }
+    case IBoundaryConditionDescriptor::Type::dirichlet: {
+      const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+        dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+      if constexpr (Dimension > 1) {
+        for (size_t i_ref_face_list = 0;
+             i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+          const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+          const RefId& ref          = ref_face_list.refId();
+          if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+            const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+            MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
+
+            Array<const FaceId> face_list = ref_face_list.list();
+            Array<const double> value_list =
+              InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+                                                                                                        mesh_data.xl(),
+                                                                                                        face_list);
+            boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
+          }
+        }
+      }
+      break;
+    }
+    case IBoundaryConditionDescriptor::Type::fourier: {
+      throw NotImplementedError("NIY");
+      break;
+    }
+    case IBoundaryConditionDescriptor::Type::neumann: {
+      const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
+        dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
+
+      for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
+           ++i_ref_face_list) {
+        const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+        const RefId& ref          = ref_face_list.refId();
+        if (ref == neumann_bc_descriptor.boundaryDescriptor()) {
+          const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension > 1) {
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+            Array<const FaceId> face_list = ref_face_list.list();
+            Array<const double> value_list =
+              InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+                                                                                                        mesh_data.xl(),
+                                                                                                        face_list);
+            boundary_condition_list.push_back(NeumannBoundaryCondition{face_list, value_list});
+          }
+        }
+      }
+      break;
+    }
+    default: {
+      is_valid_boundary_condition = false;
+    }
+    }
+    if (not is_valid_boundary_condition) {
+      std::ostringstream error_msg;
+      error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation";
+      throw NormalError(error_msg.str());
+    }
+  }
+
+  if constexpr (Dimension > 1) {
+    const CellValue<const size_t> cell_dof_number = [&] {
+      CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
+      return compute_cell_dof_number;
+    }();
+    size_t number_of_dof = mesh->numberOfCells();
+
+    const FaceValue<const size_t> face_dof_number = [&] {
+      FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+      compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+      for (const auto& boundary_condition : boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using T = std::decay_t<decltype(bc)>;
+            if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
+                          (std::is_same_v<T, FourierBoundaryCondition>) or
+                          (std::is_same_v<T, SymmetryBoundaryCondition>) or
+                          (std::is_same_v<T, DirichletBoundaryCondition>)) {
+              const auto& face_list = bc.faceList();
+
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                const FaceId face_id = face_list[i_face];
+                if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+                  std::ostringstream os;
+                  os << "The face " << face_id << " is used at least twice for boundary conditions";
+                  throw NormalError(os.str());
+                } else {
+                  compute_face_dof_number[face_id] = number_of_dof++;
+                }
+              }
+            }
+          },
+          boundary_condition);
+      }
+
+      return compute_face_dof_number;
+    }();
+
+    const FaceValue<const bool> primal_face_is_neumann = [&] {
+      FaceValue<bool> face_is_neumann{mesh->connectivity()};
+      face_is_neumann.fill(false);
+      for (const auto& boundary_condition : boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using T = std::decay_t<decltype(bc)>;
+            if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
+                          (std::is_same_v<T, FourierBoundaryCondition>) or
+                          (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+              const auto& face_list = bc.faceList();
+
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                const FaceId face_id     = face_list[i_face];
+                face_is_neumann[face_id] = true;
+              }
+            }
+          },
+          boundary_condition);
+      }
+
+      return face_is_neumann;
+    }();
+
+    const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+    const auto& face_to_cell_matrix        = mesh->connectivity().faceToCellMatrix();
+    NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
+    if (parallel::size() > 1) {
+      throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+    }
+
+    primal_node_is_on_boundary.fill(false);
+    for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+      if (face_to_cell_matrix[face_id].size() == 1) {
+        for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+          NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
+          primal_node_is_on_boundary[node_id] = true;
+        }
+      }
+    }
+
+    FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+    if (parallel::size() > 1) {
+      throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+    }
+
+    primal_face_is_on_boundary.fill(false);
+    for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+      if (face_to_cell_matrix[face_id].size() == 1) {
+        primal_face_is_on_boundary[face_id] = true;
+      }
+    }
+
+    FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
+    if (parallel::size() > 1) {
+      throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+    }
+
+    primal_face_is_dirichlet.fill(false);
+    for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+      primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
+    }
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    CellValue<double> Tj =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
+    CellValue<double> Temperature{mesh->connectivity()};
+    Temperature.fill(1);
+    FaceValue<double> Tl =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
+    FaceValue<double> Temperature_face =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
+    NodeValue<double> Tr(mesh->connectivity());
+    {
+      VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
+
+      vtk_writer.write(mesh,
+                       {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
+                        NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                       0, true);   // forces last output
+    }
+    ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face);
+    {
+      Vector<double> error{mesh->numberOfCells()};
+      CellValue<double> cell_error{mesh->connectivity()};
+      Vector<double> face_error{mesh->numberOfFaces()};
+      double error_max                        = 0.;
+      size_t cell_max                         = 0;
+      const CellValue<const double> primal_Vj = mesh_data.Vj();
+      const FaceValue<const double> mes_l     = [&] {
+        if constexpr (Dimension == 1) {
+          FaceValue<double> compute_mes_l{mesh->connectivity()};
+          compute_mes_l.fill(1);
+          return compute_mes_l;
+        } else {
+          return mesh_data.ll();
+        }
+      }();
+
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+          error[cell_id]      = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+          cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
+        });
+      parallel_for(
+        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+          if (primal_face_is_on_boundary[face_id]) {
+            face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
+          } else {
+            face_error[face_id] = 0;
+          }
+        });
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
+        if (error_max < std::abs(cell_error[cell_id])) {
+          error_max = std::abs(cell_error[cell_id]);
+          cell_max  = cell_id;
+        }
+      }
+
+      std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
+      std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n";
+      std::cout << "||Error||_2 (face)= " << std::sqrt((face_error, face_error)) << "\n";
+      std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((face_error, face_error)) << "\n";
+
+      {
+        VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
+
+        vtk_writer.write(mesh,
+                         {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
+                          NamedItemValue{"Texact", Tj}, NamedItemValue{"error", cell_error}},
+                         0,
+                         true);   // forces last output
+      }
+    }
+  } else {
+    throw NotImplementedError("not done in 1d");
+  }
+}
+
+template <size_t Dimension>
+class ParabolicHeatScheme<Dimension>::DirichletBoundaryCondition
+{
+ private:
+  const Array<const double> m_value_list;
+  const Array<const FaceId> m_face_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+    : m_value_list{value_list}, m_face_list{face_list}
+  {
+    Assert(m_value_list.size() == m_face_list.size());
+  }
+
+  ~DirichletBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class ParabolicHeatScheme<Dimension>::NeumannBoundaryCondition
+{
+ private:
+  const Array<const double> m_value_list;
+  const Array<const FaceId> m_face_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+    : m_value_list{value_list}, m_face_list{face_list}
+  {
+    Assert(m_value_list.size() == m_face_list.size());
+  }
+
+  ~NeumannBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class ParabolicHeatScheme<Dimension>::FourierBoundaryCondition
+{
+ private:
+  const Array<const double> m_coef_list;
+  const Array<const double> m_value_list;
+  const Array<const FaceId> m_face_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  const Array<const double>&
+  coefList() const
+  {
+    return m_coef_list;
+  }
+
+ public:
+  FourierBoundaryCondition(const Array<const FaceId>& face_list,
+                           const Array<const double>& coef_list,
+                           const Array<const double>& value_list)
+    : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
+  {
+    Assert(m_coef_list.size() == m_face_list.size());
+    Assert(m_value_list.size() == m_face_list.size());
+  }
+
+  ~FourierBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class ParabolicHeatScheme<Dimension>::SymmetryBoundaryCondition
+{
+ private:
+  const Array<const FaceId> m_face_list;
+
+ public:
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+ public:
+  SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+template ParabolicHeatScheme<1>::ParabolicHeatScheme(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const double&,
+  const double&);
+
+template ParabolicHeatScheme<2>::ParabolicHeatScheme(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const double&,
+  const double&);
+
+template ParabolicHeatScheme<3>::ParabolicHeatScheme(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const double&,
+  const double&);

src/language/algorithms/ParabolicHeat.hpp

+#ifndef PARABOLIC_HEAT_HPP
+#define PARABOLIC_HEAT_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+class ParabolicHeatScheme
+{
+ private:
+  class DirichletBoundaryCondition;
+  class FourierBoundaryCondition;
+  class NeumannBoundaryCondition;
+  class SymmetryBoundaryCondition;
+
+ public:
+  ParabolicHeatScheme(std::shared_ptr<const IMesh> i_mesh,
+                      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                      const FunctionSymbolId& T_id,
+                      const FunctionSymbolId& kappa_id,
+                      const FunctionSymbolId& f_id,
+                      const double& final_time,
+                      const double& dt);
+};
+
+#endif   // HEAT_DIAMOND_ALGORITHM_HPP

src/language/modules/SchemeModule.cpp

@@ -4,6 +4,8 @@
 #include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
 #include <language/algorithms/Heat5PointsAlgorithm.hpp>
 #include <language/algorithms/HeatDiamondAlgorithm.hpp>
+#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
+#include <language/algorithms/ParabolicHeat.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Mesh.hpp>
@@ -260,6 +262,74 @@ SchemeModule::SchemeModule()
 
                                       ));
 
+  this->_addBuiltinFunction("parabolicheat",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IMesh>,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
+                                   const double&, const double&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
+                                 const FunctionSymbolId& f_id, const double& final_time, const double& dt) -> void {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+                                                         f_id,   final_time,         dt};
+                                  break;
+                                }
+                                case 2: {
+                                  ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+                                                         f_id,   final_time,         dt};
+                                  break;
+                                }
+                                case 3: {
+                                  ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+                                                         f_id,   final_time,         dt};
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("heat2",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IMesh>,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
+                                 const FunctionSymbolId& f_id) -> void {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  HeatDiamondScheme2<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+                                  break;
+                                }
+                                case 2: {
+                                  HeatDiamondScheme2<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+                                  break;
+                                }
+                                case 3: {
+                                  HeatDiamondScheme2<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("elasticity",
                             std::make_shared<BuiltinFunctionEmbedder<
                               void(std::shared_ptr<const IMesh>,