"src/language/node_processor/BuiltinFunctionProcessor.hpp" did not exist on "3765c811be5acb62cda7daace60f61019ac99577"
Select Git revision
SchemeModule.cpp
-
Stéphane Del Pino authored
The aim of this singleton is to store (and eventually compute) various data associated to a given mesh: volumes, corner vectors, centroids, isobarycenters. These data will be stored as long as the mesh lives and will be computed on demand This will also allow to define really constant meshes (removing the remaining mutable node coordinates ...)
Stéphane Del Pino authoredThe aim of this singleton is to store (and eventually compute) various data associated to a given mesh: volumes, corner vectors, centroids, isobarycenters. These data will be stored as long as the mesh lives and will be computed on demand This will also allow to define really constant meshes (removing the remaining mutable node coordinates ...)
SchemeModule.cpp 11.42 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 = MeshDataLegacy<MeshType>;
using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& m_mesh;
GlaceScheme(const IMesh& 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{dynamic_cast<const MeshType&>(mesh)}
{
MeshDataType mesh_data(m_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(mesh_data);
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<MeshDataType> acoustic_solver(mesh_data, bc_list);
const CellValue<const double>& Vj = mesh_data.Vj();
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();
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)) {
vtk_writer.write(m_mesh,
{NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
NamedItemValue{"coords", m_mesh.xr()},
NamedItemValue{"cell_owner", m_mesh.connectivity().cellOwner()},
NamedItemValue{"node_owner", m_mesh.connectivity().nodeOwner()}},
t);
double dt = 0.4 * acoustic_solver.acoustic_dt(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';
acoustic_solver.computeNextStep(t, 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';
vtk_writer.write(m_mesh,
{NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
NamedItemValue{"coords", m_mesh.xr()},
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");
}
}
}
));
}