Skip to content
Snippets Groups Projects
Commit 24e89d4a authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add 3 FunctionId parameters to the glace function

These functions set initial density, velocity and pressure
parent eec76edb
No related branches found
No related tags found
1 merge request!37Feature/language
......@@ -15,9 +15,6 @@
#include <Kokkos_Core.hpp>
#include <array>
#include <cstdio>
template <typename T>
class MeshTransformation;
template <typename OutputType, typename... InputType>
......
......@@ -14,8 +14,52 @@
/////////// 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)
{
const auto flatten_args = Adapter::getFlattenArgs(function_symbol_id);
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, &flatten_args, &tokens](ItemId i) {
const int32_t t = tokens.acquire();
auto& execution_policy = context_list[t];
Adapter::convertArgs(execution_policy.currentContext(), flatten_args, 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
{
......@@ -27,7 +71,10 @@ struct GlaceScheme
const MeshType& m_mesh;
GlaceScheme(const IMesh& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
FunctionSymbolId rho_id,
FunctionSymbolId u_id,
FunctionSymbolId p_id)
: m_mesh{dynamic_cast<const MeshType&>(mesh)}
{
MeshDataType mesh_data(m_mesh);
......@@ -71,7 +118,18 @@ struct GlaceScheme
}
UnknownsType unknowns(mesh_data);
unknowns.initializeSod();
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);
......@@ -88,8 +146,22 @@ struct GlaceScheme
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);
......@@ -170,28 +242,30 @@ SchemeModule::SchemeModule()
));
this
->_addBuiltinFunction("glace",
this->_addBuiltinFunction("glace",
std::make_shared<
BuiltinFunctionEmbedder<void, std::shared_ptr<const IMesh>,
std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>>>(
std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>,
FunctionSymbolId, FunctionSymbolId, FunctionSymbolId>>(
std::function<void(std::shared_ptr<const IMesh>,
std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>)>{
std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>,
FunctionSymbolId, FunctionSymbolId, FunctionSymbolId)>{
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> void {
bc_descriptor_list,
FunctionSymbolId rho_id, FunctionSymbolId u_id, FunctionSymbolId p_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
GlaceScheme<1>{*p_mesh, bc_descriptor_list};
GlaceScheme<1>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
break;
}
case 2: {
GlaceScheme<2>{*p_mesh, bc_descriptor_list};
GlaceScheme<2>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
break;
}
case 3: {
GlaceScheme<3>{*p_mesh, bc_descriptor_list};
GlaceScheme<3>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
break;
}
default: {
......
......@@ -198,13 +198,24 @@ class PugsFunctionAdapter<OutputType(InputType...)>
[[nodiscard]] PUGS_INLINE static std::function<OutputType(DataVariant&& result)>
getResultConverter(const ASTNodeDataType& data_type)
{
if constexpr (is_tiny_vector_v<OutputType>) {
switch (data_type) {
case ASTNodeDataType::list_t: {
return [](DataVariant&& result) -> OutputType {
AggregateDataVariant& v = std::get<AggregateDataVariant>(result);
OutputType x;
for (size_t i = 0; i < x.dimension(); ++i) {
x[i] = std::get<double>(v[i]);
std::visit(
[&](auto&& vi) {
using Vi_T = std::decay_t<decltype(vi)>;
if constexpr (std::is_arithmetic_v<Vi_T>) {
x[i] = vi;
} else {
throw UnexpectedError("expecting arithmetic value");
}
},
v[i]);
}
return x;
};
......@@ -212,6 +223,33 @@ class PugsFunctionAdapter<OutputType(InputType...)>
case ASTNodeDataType::vector_t: {
return [](DataVariant&& result) -> OutputType { return std::get<OutputType>(result); };
}
case ASTNodeDataType::bool_t: {
if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
return
[](DataVariant&& result) -> OutputType { return OutputType{static_cast<double>(std::get<bool>(result))}; };
} else {
throw UnexpectedError("unexpected data_type");
}
}
case ASTNodeDataType::unsigned_int_t: {
if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
return [](DataVariant&& result) -> OutputType {
return OutputType(static_cast<double>(std::get<uint64_t>(result)));
};
} else {
throw UnexpectedError("unexpected data_type");
}
}
case ASTNodeDataType::int_t: {
if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
return [](DataVariant&& result) -> OutputType {
return OutputType{static_cast<double>(std::get<int64_t>(result))};
};
} else {
// If this point is reached must be a 0 vector
return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
}
}
case ASTNodeDataType::double_t: {
if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
return [](DataVariant&& result) -> OutputType { return OutputType{std::get<double>(result)}; };
......@@ -223,6 +261,27 @@ class PugsFunctionAdapter<OutputType(InputType...)>
throw UnexpectedError("unexpected data_type");
}
}
} else if constexpr (std::is_arithmetic_v<OutputType>) {
switch (data_type) {
case ASTNodeDataType::bool_t: {
return [](DataVariant&& result) -> OutputType { return std::get<bool>(result); };
}
case ASTNodeDataType::unsigned_int_t: {
return [](DataVariant&& result) -> OutputType { return std::get<uint64_t>(result); };
}
case ASTNodeDataType::int_t: {
return [](DataVariant&& result) -> OutputType { return std::get<int64_t>(result); };
}
case ASTNodeDataType::double_t: {
return [](DataVariant&& result) -> OutputType { return std::get<double>(result); };
}
default: {
throw UnexpectedError("unexpected data_type");
}
}
} else {
static_assert(std::is_arithmetic_v<OutputType>, "unexpected output type");
}
}
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment