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

Add local wavespeeds kinetic lagrangian scheme for multiD problems

parent f65ffd61
No related branches found
No related tags found
No related merge requests found
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include <scheme/EulerKineticSolverFirstOrderFV.hpp> #include <scheme/EulerKineticSolverFirstOrderFV.hpp>
#include <scheme/EulerKineticSolverLagrangeFV.hpp> #include <scheme/EulerKineticSolverLagrangeFV.hpp>
#include <scheme/EulerKineticSolverLagrangeMultiD.hpp> #include <scheme/EulerKineticSolverLagrangeMultiD.hpp>
#include <scheme/EulerKineticSolverLagrangeMultiDLocal.hpp>
#include <scheme/EulerKineticSolverMeanFluxMood.hpp> #include <scheme/EulerKineticSolverMeanFluxMood.hpp>
#include <scheme/EulerKineticSolverMoodFD.hpp> #include <scheme/EulerKineticSolverMoodFD.hpp>
#include <scheme/EulerKineticSolverMoodFV.hpp> #include <scheme/EulerKineticSolverMoodFV.hpp>
...@@ -860,6 +861,27 @@ KineticSchemeModule::KineticSchemeModule() ...@@ -860,6 +861,27 @@ KineticSchemeModule::KineticSchemeModule()
)); ));
this->_addBuiltinFunction("euler_kinetic_lagrange_multiD_local",
std::function(
[](const double& dt, const size_t& L, const size_t& M, const size_t& k,
const double& gamma, const double& eps, const size_t& space_order,
const size_t& time_order, const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return eulerKineticSolverLagrangeMultiDLocal(dt, L, M, k, gamma, eps, space_order,
time_order, rho, u, E, p,
bc_descriptor_list);
}
));
this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV", this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
std::function( std::function(
......
...@@ -18,6 +18,7 @@ add_library( ...@@ -18,6 +18,7 @@ add_library(
WavesEquationNonUniformCellKineticSolver.cpp WavesEquationNonUniformCellKineticSolver.cpp
EulerKineticSolver.cpp EulerKineticSolver.cpp
EulerKineticSolverLagrangeMultiD.cpp EulerKineticSolverLagrangeMultiD.cpp
EulerKineticSolverLagrangeMultiDLocal.cpp
EulerKineticSolverMeanFluxMood.cpp EulerKineticSolverMeanFluxMood.cpp
EulerKineticSolverOneFluxMood.cpp EulerKineticSolverOneFluxMood.cpp
EulerKineticSolverMoodFD.cpp EulerKineticSolverMoodFD.cpp
......
...@@ -185,7 +185,7 @@ class EulerKineticSolverLagrangeMultiD ...@@ -185,7 +185,7 @@ class EulerKineticSolverLagrangeMultiD
const size_t nb_waves = m_lambda_vector.size(); const size_t nb_waves = m_lambda_vector.size();
CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves}; CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves};
TinyVector<Dimension> inv_S = zero; TinyVector<Dimension> inv_S = zero;
TinyMatrix<Dimension> I = identity; // TinyMatrix<Dimension> I = identity;
for (size_t d = 0; d < Dimension; ++d) { for (size_t d = 0; d < Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) { for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d]; inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
...@@ -197,10 +197,8 @@ class EulerKineticSolverLagrangeMultiD ...@@ -197,10 +197,8 @@ class EulerKineticSolverLagrangeMultiD
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
M[cell_id][i_wave] = 1. / nb_waves * u[cell_id]; M[cell_id][i_wave] = 1. / nb_waves * u[cell_id];
for (size_t d1 = 0; d1 < Dimension; ++d1) { for (size_t d = 0; d < Dimension; ++d) {
for (size_t d2 = 0; d2 < Dimension; ++d2) { M[cell_id][i_wave][d] += 1. / Dimension * inv_S[d] * m_lambda_vector[i_wave][d] * p[cell_id];
M[cell_id][i_wave][d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * p[cell_id] * I(d2, d1);
}
} }
} }
}); });
...@@ -251,6 +249,158 @@ class EulerKineticSolverLagrangeMultiD ...@@ -251,6 +249,158 @@ class EulerKineticSolverLagrangeMultiD
return Fr; return Fr;
} }
NodeArray<double>
compute_Flux_Fp_lambdaj(const CellValue<const double>& p,
const CellValue<const TinyVector<Dimension>>& u,
const CellValue<const double>& lambda_j,
const NodeValue<const bool>& is_boundary_node)
{
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const size_t nb_waves = m_lambda_vector.size();
NodeArray<double> Fr(m_mesh.connectivity(), nb_waves);
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
Fr.fill(0.);
if (Dimension > 1) {
throw NotImplementedError("Non local lambda not implemented for dimension > 1");
}
// std::cout << lambda_j << "\n";
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (not is_boundary_node[node_id]) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
double lambda_r = 0.;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
lambda_r = std::max(lambda_r, lambda_j[cell_id]);
}
// Only in 1D
std::vector<TinyVector<Dimension>> lambda_r_vector(2);
lambda_r_vector[0][0] = lambda_r;
lambda_r_vector[1][0] = -lambda_r;
// std::cout << lambda_r_vector << "\n";
TinyVector<Dimension> inv_S = zero;
for (size_t d = 0; d < Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += lambda_r_vector[i][d] * lambda_r_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
double sum_li_njr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
double li_njr = dot(lambda_r_vector[i_wave], njr(cell_id, i_node));
if (li_njr > 1e-14) {
double Fj = (1. / nb_waves) * p[cell_id];
for (size_t d = 0; d < Dimension; ++d) {
Fj += inv_S[d] * lambda_r_vector[i_wave][d] * (lambda_r * lambda_r * u[cell_id][d]);
}
Fr[node_id][i_wave] += li_njr * Fj;
sum_li_njr += li_njr;
}
}
if (sum_li_njr != 0) {
Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
}
}
}
});
return Fr;
}
NodeArray<TinyVector<Dimension>>
compute_Flux_Fu_lambdaj(const CellValue<const double>& p,
const CellValue<const TinyVector<Dimension>>& u,
const CellValue<const double>& lambda_j,
const NodeValue<const bool>& is_boundary_node)
{
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const size_t nb_waves = m_lambda_vector.size();
NodeArray<TinyVector<Dimension>> Fr(m_mesh.connectivity(), nb_waves);
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
Fr.fill(zero);
if (Dimension > 1) {
throw NotImplementedError("Non local lambda not implemented for dimension > 1");
}
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (not is_boundary_node[node_id]) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
double lambda_r = 0.;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
lambda_r = std::max(lambda_r, lambda_j[cell_id]);
}
// Only in 1D
std::vector<TinyVector<Dimension>> lambda_r_vector(2);
lambda_r_vector[0][0] = lambda_r;
lambda_r_vector[1][0] = -lambda_r;
TinyVector<Dimension> inv_S = zero;
for (size_t d = 0; d < Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += lambda_r_vector[i][d] * lambda_r_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
double sum_li_njr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
double li_njr = dot(lambda_r_vector[i_wave], njr(cell_id, i_node));
if (li_njr > 1e-14) {
TinyVector<Dimension> Fj = (1. / nb_waves) * u[cell_id];
for (size_t d1 = 0; d1 < Dimension; ++d1) {
for (size_t d2 = 0; d2 < Dimension; ++d2) {
TinyMatrix<Dimension> I = identity;
Fj[d1] += inv_S[d2] * lambda_r_vector[i_wave][d2] * p[cell_id] * I(d2, d1);
}
}
Fr[node_id][i_wave] += li_njr * Fj;
sum_li_njr += li_njr;
}
}
if (sum_li_njr != 0) {
Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
}
}
}
});
return Fr;
}
void void
apply_bc(NodeArray<double>& Fr_p, apply_bc(NodeArray<double>& Fr_p,
NodeArray<TinyVector<Dimension>>& Fr_u, NodeArray<TinyVector<Dimension>>& Fr_u,
...@@ -503,6 +653,39 @@ class EulerKineticSolverLagrangeMultiD ...@@ -503,6 +653,39 @@ class EulerKineticSolverLagrangeMultiD
return u; return u;
} }
const NodeValue<TinyVector<Dimension>>
compute_velocity_flux_lambdaj(NodeArray<double>& F, const CellValue<const double>& lambda_j)
{
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
NodeValue<TinyVector<Dimension>> u{p_mesh->connectivity()};
const size_t nb_waves = m_lambda_vector.size();
u.fill(zero);
parallel_for(
m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
double lambda_r = 0.;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
lambda_r = std::max(lambda_r, lambda_j[cell_id]);
}
const double inv_lambda_squared = (1. / (lambda_r * lambda_r));
std::vector<TinyVector<Dimension>> lambda_r_vector(2);
lambda_r_vector[0][0] = lambda_r;
lambda_r_vector[1][0] = -lambda_r;
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
u[node_id] += F[node_id][i_wave] * lambda_r_vector[i_wave];
}
u[node_id] = inv_lambda_squared * u[node_id];
});
return u;
}
const NodeValue<TinyVector<Dimension>> const NodeValue<TinyVector<Dimension>>
compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F) compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F)
{ {
...@@ -514,7 +697,41 @@ class EulerKineticSolverLagrangeMultiD ...@@ -514,7 +697,41 @@ class EulerKineticSolverLagrangeMultiD
m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t dim = 0; dim < Dimension; ++dim) { for (size_t dim = 0; dim < Dimension; ++dim) {
p[node_id][dim] += F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim]; p[node_id][dim] += dot(F[node_id][i_wave], m_lambda_vector[i_wave]);
// p[node_id][0] += 0.5 * F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
// p[node_id][1] += 0.5 * F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
}
}
});
return p;
}
const NodeValue<TinyVector<Dimension>>
compute_pressure_flux_lambdaj(NodeArray<TinyVector<Dimension>>& F, const CellValue<const double>& lambda_j)
{
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
NodeValue<TinyVector<Dimension>> p{p_mesh->connectivity()};
const size_t nb_waves = m_lambda_vector.size();
p.fill(zero);
parallel_for(
m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
double lambda_r = 0.;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
lambda_r = std::max(lambda_r, lambda_j[cell_id]);
}
std::vector<TinyVector<Dimension>> lambda_r_vector(2);
lambda_r_vector[0][0] = lambda_r;
lambda_r_vector[1][0] = -lambda_r;
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t dim = 0; dim < Dimension; ++dim) {
p[node_id][dim] += F[node_id][i_wave][dim] * lambda_r_vector[i_wave][dim];
} }
} }
}); });
...@@ -638,14 +855,37 @@ class EulerKineticSolverLagrangeMultiD ...@@ -638,14 +855,37 @@ class EulerKineticSolverLagrangeMultiD
CellValue<double> p_np1 = copy(p); CellValue<double> p_np1 = copy(p);
CellValue<double> E_np1 = copy(E); CellValue<double> E_np1 = copy(E);
const CellValue<const double> cj = [&] {
CellValue<double> computed_cj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(),
PUGS_LAMBDA(const CellId cell_id) { computed_cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]); });
return computed_cj;
}();
const CellValue<const double> lambda_j = [&] {
CellValue<double> computed_lambda_j{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(),
PUGS_LAMBDA(const CellId cell_id) { computed_lambda_j[cell_id] = std::abs(rho[cell_id] * cj[cell_id]); });
return computed_lambda_j;
}();
// std::cout << "c = " << cj << ", \n lambda = " << lambda_j << ", \n rho = " << rho << "\n";
const CellArray<const double> F_p = compute_M_p(p, u, lambda); const CellArray<const double> F_p = compute_M_p(p, u, lambda);
const CellArray<const TinyVector<Dimension>> F_u = compute_M_u(u, p); const CellArray<const TinyVector<Dimension>> F_u = compute_M_u(u, p);
NodeArray<double> Fr_p = compute_Flux<double>(F_p, is_boundary_node); NodeArray<double> Fr_p = compute_Flux<double>(F_p, is_boundary_node);
NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux<TinyVector<Dimension>>(F_u, is_boundary_node); NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux<TinyVector<Dimension>>(F_u, is_boundary_node);
// NodeArray<double> Fr_p = compute_Flux_Fp_lambdaj(p, u, lambda_j, is_boundary_node);
// NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux_Fu_lambdaj(p, u, lambda_j, is_boundary_node);
apply_bc(Fr_p, Fr_u, p, u, F_p, F_u, lambda, bc_list); apply_bc(Fr_p, Fr_u, p, u, F_p, F_u, lambda, bc_list);
// const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux_lambdaj(Fr_p, lambda_j);
// const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux_lambdaj(Fr_u, lambda_j);
const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda); const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda);
const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u); const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u);
......
This diff is collapsed.
#ifndef EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
#define EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
#include <language/utils/FunctionSymbolId.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
class IBoundaryConditionDescriptor;
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
eulerKineticSolverLagrangeMultiDLocal(
const double& dt,
const size_t& L,
const size_t& M,
const size_t& k,
const double& gamma,
const double& eps,
const size_t& space_order,
const size_t& time_order,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
#endif // EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment