Skip to content
Snippets Groups Projects
Select Git revision
  • f38904e3d21ff1b6cc44d9c7fda05ac1961a0d59
  • develop default protected
  • feature/gmsh-reader
  • 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
  • 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

SchemeModule.cpp

Blame
    • Stéphane Del Pino's avatar
      f38904e3
      Fix MeshData and MeshDataManager · f38904e3
      Stéphane Del Pino authored
      MeshData constructor can now only be called by MeshDataManager! This
      ensures only one construction of the data storage for each mesh
      
      Also fix following issues
      - volume calculation and checking was dangerously nested
      - few iterators were poorly named (due to copy/paste in MeshDataManager)
      f38904e3
      History
      Fix MeshData and MeshDataManager
      Stéphane Del Pino authored
      MeshData constructor can now only be called by MeshDataManager! This
      ensures only one construction of the data storage for each mesh
      
      Also fix following issues
      - volume calculation and checking was dangerously nested
      - few iterators were poorly named (due to copy/paste in MeshDataManager)
    SchemeModule.cpp 11.77 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/SymmetryBoundaryConditionDescriptor.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 <size_t Dimension>
    struct GlaceScheme
    {
      using ConnectivityType = Connectivity<Dimension>;
      using MeshType         = Mesh<ConnectivityType>;
      using MeshDataType     = MeshData<Dimension>;
      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
    
      std::shared_ptr<const MeshType> m_mesh;
    
      GlaceScheme(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';
    
        std::vector<BoundaryConditionHandler> 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: {
              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()) {
                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
                    new SymmetryBoundaryCondition<MeshType::Dimension>(
                      MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list));
                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
                  bc_list.push_back(BoundaryConditionHandler(bc));
                }
              }
              break;
            }
            default: {
              throw UnexpectedError("Unknown BCDescription\n");
            }
            }
          }
        }
    
        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);
    
        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("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: {
                                      GlaceScheme<1>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                      break;
                                    }
                                    case 2: {
                                      GlaceScheme<2>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                      break;
                                    }
                                    case 3: {
                                      GlaceScheme<3>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
                                      break;
                                    }
                                    default: {
                                      throw UnexpectedError("invalid mesh dimension");
                                    }
                                    }
                                  }
    
                                  ));
    }