Skip to content
Snippets Groups Projects
Select Git revision
  • f5fa77f1f0e8f6fe36c628b9512b31e4c274221d
  • develop default protected
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • feature/hypoelasticity-clean
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

OperatorRepository.hpp

Blame
  • SchemeModule.cpp 18.22 KiB
    #include <language/modules/SchemeModule.hpp>
    
    #include <language/utils/BuiltinFunctionEmbedder.hpp>
    #include <language/utils/TypeDescriptor.hpp>
    #include <mesh/Mesh.hpp>
    #include <scheme/AcousticSolver.hpp>
    #include <scheme/IBoundaryConditionDescriptor.hpp>
    #include <scheme/IBoundaryDescriptor.hpp>
    #include <scheme/NamedBoundaryDescriptor.hpp>
    #include <scheme/NumberedBoundaryDescriptor.hpp>
    #include <scheme/PressureBoundaryConditionDescriptor.hpp>
    #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
    #include <scheme/VelocityBoundaryConditionDescriptor.hpp>
    
    #include <memory>
    
    /////////// TEMPORARY
    
    #include <language/utils/PugsFunctionAdapter.hpp>
    #include <output/VTKWriter.hpp>
    
    template <typename T>
    class InterpolateItemValue;
    template <typename OutputType, typename InputType>
    class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
    {
      static constexpr size_t Dimension = OutputType::Dimension;
      using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
    
     public:
      template <ItemType item_type>
      static inline ItemValue<OutputType, item_type>
      interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
      {
        auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
        auto convert_result = Adapter::getResultConverter(expression.m_data_type);
    
        Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
    
        using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
        Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
        const IConnectivity& connectivity = *position.connectivity_ptr();
    
        ItemValue<OutputType, item_type> value(connectivity);
        using ItemId = ItemIdT<item_type>;
    
        parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
          const int32_t t = tokens.acquire();
    
          auto& execution_policy = context_list[t];
    
          Adapter::convertArgs(execution_policy.currentContext(), position[i]);
          auto result = expression.execute(execution_policy);
          value[i]    = convert_result(std::move(result));
    
          tokens.release(t);
        });
    
        return value;
      }
    
      template <ItemType item_type>
      static inline Array<OutputType>
      interpolate(const FunctionSymbolId& function_symbol_id,
                  const ItemValue<const InputType, item_type>& position,
                  const Array<const ItemIdT<item_type>>& list_of_items)
      {
        auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
        auto convert_result = Adapter::getResultConverter(expression.m_data_type);
    
        Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
    
        using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
        Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
    
        Array<OutputType> value{list_of_items.size()};
        using ItemId = ItemIdT<item_type>;
    
        parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) {
          ItemId item_id  = list_of_items[i_item];
          const int32_t t = tokens.acquire();
    
          auto& execution_policy = context_list[t];
    
          Adapter::convertArgs(execution_policy.currentContext(), position[item_id]);
          auto result   = expression.execute(execution_policy);
          value[i_item] = convert_result(std::move(result));
    
          tokens.release(t);
        });
    
        return value;
      }
    };
    
    template <size_t Dimension, AcousticSolverType solver_type>
    struct AcousticSolverScheme
    {
      using ConnectivityType = Connectivity<Dimension>;
      using MeshType         = Mesh<ConnectivityType>;
      using MeshDataType     = MeshData<Dimension>;
      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
    
      std::shared_ptr<const MeshType> m_mesh;
    
      AcousticSolverScheme(std::shared_ptr<const IMesh> i_mesh,
                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                           const FunctionSymbolId& rho_id,
                           const FunctionSymbolId& u_id,
                           const FunctionSymbolId& p_id)
        : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
      {
        std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
    
        typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
        {
          constexpr ItemType FaceType = [] {
            if constexpr (Dimension > 1) {
              return ItemType::face;
            } else {
              return ItemType::node;
            }
          }();
    
          for (const auto& bc_descriptor : bc_descriptor_list) {
            switch (bc_descriptor->type()) {
            case IBoundaryConditionDescriptor::Type::symmetry: {
              using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
    
              const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
                dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
              for (size_t i_ref_face_list = 0;
                   i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
                const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
                const RefId& ref          = ref_face_list.refId();
                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                  bc_list.push_back(
                    SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list)});
                }
              }
              break;
            }
            case IBoundaryConditionDescriptor::Type::velocity: {
              using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
    
              const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
                dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor);
              for (size_t i_ref_face_list = 0;
                   i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
                const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
                const RefId& ref          = ref_face_list.refId();
                if (ref == velocity_bc_descriptor.boundaryDescriptor()) {
                  MeshNodeBoundary<Dimension> mesh_node_boundary{m_mesh, ref_face_list};
    
                  const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
    
                  const auto& node_list = mesh_node_boundary.nodeList();
    
                  Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
                    TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
    
                  bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
                }
              }
              break;
            }
            case IBoundaryConditionDescriptor::Type::pressure: {
              using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
    
              const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
                dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor);
              for (size_t i_ref_face_list = 0;
                   i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
                const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
                const RefId& ref          = ref_face_list.refId();
                if (ref == pressure_bc_descriptor.boundaryDescriptor()) {
                  const auto& face_list = ref_face_list.list();
    
                  const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId();
    
                  Array<const double> face_values = [&] {
                    if constexpr (Dimension == 1) {
                      return InterpolateItemValue<double(
                        TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, m_mesh->xr(), face_list);
                    } else {
                      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
    
                      return InterpolateItemValue<double(
                        TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
                    }
                  }();
    
                  bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
                }
              }
              break;
            }
            default: {
              std::ostringstream error_msg;
              error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
              throw NormalError(error_msg.str());
            }
            }
          }
        }
    
        UnknownsType unknowns(*m_mesh);
    
        {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
    
          unknowns.rhoj() =
            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
                                                                                                      mesh_data.xj());
    
          unknowns.pj() =
            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
    
          unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>(
            TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj());
        }
        unknowns.gammaj().fill(1.4);
    
        AcousticSolver acoustic_solver(m_mesh, bc_list, solver_type);
    
        const double tmax = 0.2;
        double t          = 0;
    
        int itermax   = std::numeric_limits<int>::max();
        int iteration = 0;
    
        CellValue<double>& rhoj              = unknowns.rhoj();
        CellValue<double>& ej                = unknowns.ej();
        CellValue<double>& pj                = unknowns.pj();
        CellValue<double>& gammaj            = unknowns.gammaj();
        CellValue<double>& cj                = unknowns.cj();
        CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
        CellValue<double>& Ej                = unknowns.Ej();
        CellValue<double>& mj                = unknowns.mj();
        CellValue<double>& inv_mj            = unknowns.invMj();
    
        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
        block_eos.updateEandCFromRhoP();
    
        {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
    
          const CellValue<const double> Vj = mesh_data.Vj();
    
          parallel_for(
            m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
    
          parallel_for(
            m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
    
          parallel_for(
            m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
        }
    
        VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
    
        while ((t < tmax) and (iteration < itermax)) {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
    
          vtk_writer.write(m_mesh,
                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                            NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
                            NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
                            NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
                           t);
    
          double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
          if (t + dt > tmax) {
            dt = tmax - t;
          }
    
          std::cout.setf(std::cout.scientific);
          std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
                    << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
                    << rang::style::reset << '\n';
    
          m_mesh = acoustic_solver.computeNextStep(dt, unknowns);
    
          block_eos.updatePandCFromRhoE();
    
          t += dt;
          ++iteration;
        }
        std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
                  << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
                  << rang::style::reset << '\n';
        {
          MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
    
          vtk_writer.write(m_mesh,
                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                            NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
                            NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
                            NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
                           t, true);   // forces last output
        }
      }
    };
    
    SchemeModule::SchemeModule()
    {
      this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
      this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>);
    
      this->_addBuiltinFunction("boundaryName",
                                std::make_shared<
                                  BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>(const std::string&)>>(
    
                                  [](const std::string& boundary_name) -> std::shared_ptr<const IBoundaryDescriptor> {
                                    return std::make_shared<NamedBoundaryDescriptor>(boundary_name);
                                  }
    
                                  ));
    
      this->_addBuiltinFunction("boundaryTag",
                                std::make_shared<
                                  BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>(int64_t)>>(
    
                                  [](int64_t boundary_tag) -> std::shared_ptr<const IBoundaryDescriptor> {
                                    return std::make_shared<NumberedBoundaryDescriptor>(boundary_tag);
                                  }
    
                                  ));
    
      this
        ->_addBuiltinFunction("symmetry",
                              std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
                                std::shared_ptr<const IBoundaryDescriptor>)>>(
    
                                [](std::shared_ptr<const IBoundaryDescriptor> boundary)
                                  -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                  return std::make_shared<SymmetryBoundaryConditionDescriptor>(boundary);
                                }
    
                                ));
    
      this->_addBuiltinFunction("pressure",
                                std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                                  const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                      const FunctionSymbolId&)>>(
    
                                  [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                     const FunctionSymbolId& pressure_id)
                                    -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                    return std::make_shared<PressureBoundaryConditionDescriptor>(boundary, pressure_id);
                                  }
    
                                  ));
    
      this->_addBuiltinFunction("velocity",
                                std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
                                  const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
                                                                      const FunctionSymbolId&)>>(
    
                                  [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                     const FunctionSymbolId& velocity_id)
                                    -> std::shared_ptr<const IBoundaryConditionDescriptor> {
                                    return std::make_shared<VelocityBoundaryConditionDescriptor>(boundary, velocity_id);
                                  }
    
                                  ));
    
      this->_addBuiltinFunction(
        "glace",
        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& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
            switch (p_mesh->dimension()) {
            case 1: {
              AcousticSolverScheme<1, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            case 2: {
              AcousticSolverScheme<2, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            case 3: {
              AcousticSolverScheme<3, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            default: {
              throw UnexpectedError("invalid mesh dimension");
            }
            }
          }
    
          ));
    
      this->_addBuiltinFunction(
        "eucclhyd",
        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& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
            switch (p_mesh->dimension()) {
            case 1: {
              AcousticSolverScheme<1, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            case 2: {
              AcousticSolverScheme<2, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            case 3: {
              AcousticSolverScheme<3, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
              break;
            }
            default: {
              throw UnexpectedError("invalid mesh dimension");
            }
            }
          }
    
          ));
    }