Skip to content
Snippets Groups Projects
Commit af488bd0 authored by clovis schoeck's avatar clovis schoeck
Browse files

Navier-Stokes GKS solver programmed

parent f543a21c
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
src/language/modules/SchemeModule.cpp
+ 3
− 3
View file @ af488bd0
......@@ -459,11 +459,11 @@ SchemeModule::SchemeModule()
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE,
const std::shared_ptr<const DiscreteFunctionVariant>& tau, const double& delta, const double& dt)
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE, const double& viscosity, const double& delta,
const double& dt)
-> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return GKSHandler{getCommonMesh({rho, rhoU, rhoE})}.solver().gksNavier(rho, rhoU, rhoE, tau, delta, dt);
return GKSHandler{getCommonMesh({rho, rhoU, rhoE})}.solver().gksNavier(rho, rhoU, rhoE, viscosity, delta, dt);
}
));
......
src/scheme/GKS2.cpp
+ 33
− 43
View file @ af488bd0
......@@ -33,7 +33,7 @@ class GKS2
if (tau_n[cell_id] == 0)
eta[cell_id] = 0;
else
eta[cell_id] = tau_n[cell_id]; //(tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
eta[cell_id] = (tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
}
// std::cout << "eta = " << eta << std::endl;
......@@ -83,6 +83,10 @@ class GKS2
CellValue<Rd> U{p_mesh->connectivity()};
// U.fill(Rd(0));
CellValue<double> erf_term_pos(mesh.connectivity());
CellValue<double> erf_term_neg(mesh.connectivity());
CellValue<double> exp_term(mesh.connectivity());
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
U[cell_id][0] = rho_U_n[cell_id][0] / rho_n[cell_id];
......@@ -90,6 +94,14 @@ class GKS2
lambda[cell_id] = 0.5 * (1. + delta) * rho_n[cell_id] / (2 * rho_E_n[cell_id] - rho_U_2);
});
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
double U_2 = U[cell_id][0] * U[cell_id][0];
erf_term_pos[cell_id] = 1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]);
erf_term_neg[cell_id] = 1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]);
exp_term[cell_id] = std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
});
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id];
......@@ -106,35 +118,22 @@ class GKS2
continue;
if (l == 0) {
double U_2_left = U[node_cells[l]][0] * U[node_cells[l]][0];
rho_cell_left =
rho_n[node_cells[l]] * (1 + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
rho_cell_left = rho_n[node_cells[l]] * erf_term_pos[node_cells[l]];
rho_U_cell_left[0] =
rho_U_n[node_cells[l]][0] * (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
rho_n[node_cells[l]] * std::exp(-lambda[node_cells[l]] * U_2_left) /
std::sqrt(pi * lambda[node_cells[l]]);
rho_E_cell_left =
rho_E_n[node_cells[l]] * (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
0.5 * rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_left) /
std::sqrt(pi * lambda[node_cells[l]]);
} else {
double U_2_right = U[node_cells[l]][0] * U[node_cells[l]][0];
rho_U_n[node_cells[l]][0] * erf_term_pos[node_cells[l]] + rho_n[node_cells[l]] * exp_term[node_cells[l]];
rho_E_cell_left = rho_E_n[node_cells[l]] * erf_term_pos[node_cells[l]] +
0.5 * rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
} else {
rho_cell_right =
rho_n[node_cells[l]] * (1 - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
rho_n[node_cells[l]] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
rho_U_cell_right[0] =
rho_U_n[node_cells[l]][0] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
rho_n[node_cells[l]] * std::exp(-lambda[node_cells[l]] * U_2_right) /
std::sqrt(pi * lambda[node_cells[l]]);
rho_E_cell_right =
rho_E_n[node_cells[l]] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
0.5 * rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_right) /
std::sqrt(pi * lambda[node_cells[l]]);
rho_U_n[node_cells[l]][0] * erf_term_neg[node_cells[l]] - rho_n[node_cells[l]] * exp_term[node_cells[l]];
rho_E_cell_right = rho_E_n[node_cells[l]] * erf_term_neg[node_cells[l]] -
0.5 * rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
}
}
rho_node[node_id] = 0.5 * (rho_cell_left + rho_cell_right);
......@@ -159,28 +158,26 @@ class GKS2
double U_2_left = U[node_cells[l]][0] * U[node_cells[l]][0];
double rho_U_2_left = rho_U_n[node_cells[l]][0] * U[node_cells[l]][0];
F2_fn_left[0] = (rho_U_2_left + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
(1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_left) /
std::sqrt(pi * lambda[node_cells[l]]);
F2_fn_left[0] =
(rho_U_2_left + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) * erf_term_pos[node_cells[l]] +
rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
F3_fn_left[0] = 0.5 * rho_U_n[node_cells[l]][0] * (U_2_left + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
(1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
erf_term_pos[node_cells[l]] +
0.5 * rho_n[node_cells[l]] * (U_2_left + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
std::exp(-lambda[node_cells[l]] * U_2_left) / std::sqrt(pi * lambda[node_cells[l]]);
exp_term[node_cells[l]];
} else {
double U_2_right = U[node_cells[l]][0] * U[node_cells[l]][0];
double rho_U_2_right = rho_U_n[node_cells[l]][0] * U[node_cells[l]][0];
F2_fn_right[0] = (rho_U_2_right + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
(1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_right) /
std::sqrt(pi * lambda[node_cells[l]]);
F2_fn_right[0] =
(rho_U_2_right + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) * erf_term_neg[node_cells[l]] -
rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
F3_fn_right[0] = 0.5 * rho_U_n[node_cells[l]][0] * (U_2_right + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
(1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
erf_term_neg[node_cells[l]] -
0.5 * rho_n[node_cells[l]] * (U_2_right + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
std::exp(-lambda[node_cells[l]] * U_2_right) / std::sqrt(pi * lambda[node_cells[l]]);
exp_term[node_cells[l]];
}
}
......@@ -215,13 +212,6 @@ class GKS2
const size_t& nb_nodes = rho_flux_G.numberOfSubValues(cell_id);
for (size_t r = 0; r < nb_nodes; r++) {
rho[cell_id] -= dt / Vj[cell_id] * rho_flux_G(cell_id, r)[0];
// rho_U[cell_id] -= dt / Vj[cell_id] *
// ((1 - eta[cell_id]) * rho_U_flux_G(cell_id, r) + eta[cell_id] *
// (rho_U_flux_Ffn(cell_id,
// r)));
// rho_E[cell_id] -=
// dt / Vj[cell_id] *
// ((1 - eta[cell_id]) * rho_E_flux_G(cell_id, r)[0] + eta[cell_id] * (rho_E_flux_Ffn(cell_id, r)[0]));
rho_U[cell_id] -=
dt / Vj[cell_id] *
(rho_U_flux_G(cell_id, r) + eta[cell_id] * (rho_U_flux_Ffn(cell_id, r) - rho_U_flux_G(cell_id, r)));
......
src/scheme/GKSNavier.cpp
+ 461
− 236
View file @ af488bd0
This diff is collapsed.
src/scheme/GKSNavier.hpp
+ 2
− 9
View file @ af488bd0
......@@ -19,13 +19,6 @@ double gks_inv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c,
class GKSHandler
{
public:
enum class SolverType
{
Glace,
Eucclhyd
};
private:
struct IGKSNAVIER
{
......@@ -35,7 +28,7 @@ class GKSHandler
compute_fluxes(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
const double& viscosity,
const double& delta,
const double& dt) const = 0;
......@@ -56,7 +49,7 @@ class GKSHandler
gksNavier(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
const double& viscosity,
const double& delta,
const double& dt) const = 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment