Skip to content
Snippets Groups Projects
Commit ef580441 authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Create kinetic scheme module

parent 06149ba1
No related branches found
No related tags found
No related merge requests found
......@@ -6,6 +6,7 @@ add_library(
BuiltinModule.cpp
CoreModule.cpp
DevUtilsModule.cpp
KineticSchemeModule.cpp
LinearSolverModule.cpp
MathFunctionRegisterForVh.cpp
MathModule.cpp
......
This diff is collapsed.
#ifndef KINETIC_SCHEME_MODULE_HPP
#define KINETIC_SCHEME_MODULE_HPP
#include <language/modules/BuiltinModule.hpp>
class KineticSchemeModule : public BuiltinModule
{
public:
std::string_view
name() const final
{
return "kinetic_scheme";
}
void registerOperators() const final;
void registerCheckpointResume() const final;
KineticSchemeModule();
~KineticSchemeModule() = default;
};
#endif // KINETIC_SCHEME_MODULE_HPP
This diff is collapsed.
......@@ -1077,6 +1077,104 @@ class EulerKineticSolverMultiD
NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
nr.fill(zero);
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
nr[node_id] += njr(cell_id, i_node_cell);
}
nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
auto nxn = tensorProduct(nr[node_id], nr[node_id]);
auto txt = I - nxn;
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& node_to_face = node_to_face_matrix[node_id];
const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
double sum = 0;
Fr_rho[node_id][i_wave] = 0;
Fr_rho_u[node_id][i_wave] = zero;
Fr_rho_E[node_id][i_wave] = 0;
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
const FaceId face_id = node_to_face[i_face];
const auto& face_to_cell = face_to_cell_matrix[face_id];
const size_t i_node_face = node_local_number_in_its_face[i_face];
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
if (li_nlr > 0) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
sum += li_nlr;
}
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
if (li_nlr_prime > 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
double rho_Ej_prime = rho_E[cell_id];
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyMatrix<Dimension> A_rho_u_prime =
1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves;
TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
double Fn_rho_E_prime = rho_Ej_prime / nb_waves;
for (size_t d1 = 0; d1 < Dimension; ++d1) {
Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
for (size_t d2 = 0; d2 < Dimension; ++d2) {
Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
}
Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
}
Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime;
Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime;
Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime;
sum += li_nlr_prime;
}
}
}
if (sum != 0) {
Fr_rho[node_id][i_wave] /= sum;
Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
Fr_rho_E[node_id][i_wave] /= sum;
}
}
}
} else if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto node_list = bc.nodeList();
TinyMatrix<Dimension> I = identity;
NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
nr.fill(zero);
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment