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

GKS Navier - inverse of M with Gram-Schmidt corrected, and optimisation

parent af488bd0
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/scheme/GKSNavier.cpp
+ 66
− 126
View file @ fc8dffa8
...@@ -31,10 +31,8 @@ gks_inv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v, ...@@ -31,10 +31,8 @@ gks_inv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
CellValue<double> local_inv_dt{mesh.connectivity()}; CellValue<double> local_inv_dt{mesh.connectivity()};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { mesh.numberOfCells(),
// local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id]); PUGS_LAMBDA(CellId cell_id) { local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id]); });
local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id] * Vj[cell_id]);
});
return max(local_inv_dt); return max(local_inv_dt);
} else { } else {
...@@ -75,9 +73,9 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -75,9 +73,9 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
double p = rho / (2 * lambda); double p = rho / (2 * lambda);
M(0, 0) = rho; M(0, 0) = rho;
M(1, 1) = rhoU[0] * U[0] + p; M(1, 1) = rho * U_2 + p;
M(2, 2) = M(2, 2) =
0.25 * rho * U_2 * U_2 + (delta + 4) * 0.5 * U_2 * p + 0.125 * (delta * delta + 6 * delta + 8) * p / lambda; 0.25 * rho * U_2 * U_2 + 0.5 * (delta + 3) * U_2 * p + 0.125 * (delta * delta + 4 * delta + 3) * p / lambda;
M(1, 0) = rhoU[0]; M(1, 0) = rhoU[0];
M(0, 1) = M(1, 0); M(0, 1) = M(1, 0);
...@@ -85,11 +83,10 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -85,11 +83,10 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
M(2, 0) = rhoE; M(2, 0) = rhoE;
M(0, 2) = M(2, 0); M(0, 2) = M(2, 0);
M(2, 1) = 0.5 * rhoU[0] * U_2 + (delta + 4) * 0.5 * U[0] * p; M(2, 1) = 0.5 * (rhoU[0] * U_2 + (delta + 3) * U[0] * p);
M(1, 2) = M(2, 1); M(1, 2) = M(2, 1);
const TinyMatrix<SIZEproblem> M_bis = M; const TinyMatrix<SIZEproblem> inverseM = inverse(M);
const TinyMatrix<SIZEproblem> inverseM = inverse(M_bis);
return inverseM; return inverseM;
} }
...@@ -104,12 +101,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -104,12 +101,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
TinyMatrix<SIZEproblem> invMatrixM(zero); TinyMatrix<SIZEproblem> invMatrixM(zero);
TinyMatrix<SIZEproblem> MatrixP(zero); TinyMatrix<SIZEproblem> MatrixP(zero);
TinyMatrix<SIZEproblem> invMatrixP(zero);
TinyMatrix<SIZEproblem> invMatrixMtilde(zero); TinyMatrix<SIZEproblem> invMatrixMtilde(zero);
invMatrixMtilde(0, 0) = 1; invMatrixMtilde(0, 0) = 1;
invMatrixMtilde(1, 1) = 1; invMatrixMtilde(1, 1) = 1;
invMatrixMtilde(2, 2) = 2 * sqrt2 * lambda / (std::sqrt(delta + 1) * sqrt_rho); invMatrixMtilde(2, 2) = (8 * lambda * lambda) / (rho * (delta + 1));
MatrixP(0, 0) = 1 / sqrt_rho; MatrixP(0, 0) = 1 / sqrt_rho;
MatrixP(1, 1) = sqrt2 * sqrt_lambda / sqrt_rho; MatrixP(1, 1) = sqrt2 * sqrt_lambda / sqrt_rho;
...@@ -118,16 +114,22 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -118,16 +114,22 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
MatrixP(2, 0) = 0.5 * U_2 - 0.25 * (delta + 1) / lambda; MatrixP(2, 0) = 0.5 * U_2 - 0.25 * (delta + 1) / lambda;
MatrixP(2, 1) = -U[0]; MatrixP(2, 1) = -U[0];
invMatrixP(0, 0) = sqrt_rho; // TinyMatrix<SIZEproblem> MatrixM(zero);
invMatrixP(1, 1) = sqrt_rho / (sqrt2 * sqrt_lambda); // TinyMatrix<SIZEproblem> invMatrixP(zero);
invMatrixP(2, 2) = 1; // TinyMatrix<SIZEproblem> MatrixMtilde(zero);
invMatrixP(1, 0) = sqrt_rho * U[0];
invMatrixP(2, 0) = sqrt_rho * (0.5 * U_2 + 0.25 * (delta + 1) / lambda); // MatrixMtilde(0, 0) = 1;
invMatrixP(2, 1) = U[0] * sqrt_rho / (sqrt2 * sqrt_lambda); // MatrixMtilde(1, 1) = 1;
// MatrixMtilde(2, 2) = (rho * (delta + 1)) / (8 * lambda * lambda);
invMatrixM = MatrixP * invMatrixMtilde * invMatrixP; // invMatrixP(0, 0) = sqrt_rho;
// invMatrixP(1, 1) = sqrt_rho / (sqrt2 * sqrt_lambda);
// invMatrixP(2, 2) = 1;
// invMatrixP(1, 0) = sqrt_rho * U[0];
// invMatrixP(2, 0) = sqrt_rho * (0.5 * U_2 + 0.25 * (delta + 1) / lambda);
// invMatrixP(2, 1) = U[0] * sqrt_rho / (sqrt2 * sqrt_lambda);
// std::cout << std ::endl << invMatrixP << std ::endl << std ::endl; invMatrixM = transpose(MatrixP) * invMatrixMtilde * (MatrixP);
return invMatrixM; return invMatrixM;
} }
...@@ -160,7 +162,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -160,7 +162,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
return xi2_moments; return xi2_moments;
} }
std::tuple<CellValue<const TinyVector<SIZEproblem>>, CellValue<const TinyVector<SIZEproblem>>> CellValue<const TinyVector<SIZEproblem>>
_computeGradConservVariables_for_a(const MeshType& mesh, _computeGradConservVariables_for_a(const MeshType& mesh,
const DiscreteScalarFunction& rho, const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& rhoU, const DiscreteVectorFunction& rhoU,
...@@ -174,8 +176,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -174,8 +176,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()}; CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()};
Gradj.fill(zero); Gradj.fill(zero);
CellValue<TinyVector<SIZEproblem>> Gradjpu{mesh.connectivity()};
Gradjpu.fill(zero);
NodeValue<double> rho_node_mean{mesh.connectivity()}; NodeValue<double> rho_node_mean{mesh.connectivity()};
NodeValue<Rd> rhoU_node_mean{mesh.connectivity()}; NodeValue<Rd> rhoU_node_mean{mesh.connectivity()};
...@@ -201,23 +201,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -201,23 +201,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoU_node_mean[node_id] *= 0.5; rhoU_node_mean[node_id] *= 0.5;
rhoE_node_mean[node_id] *= 0.5; rhoE_node_mean[node_id] *= 0.5;
for (size_t l = 0; l < node_cells.size(); l++) { Gradj[node_cells[0]][0] = (rho_node_mean[node_id] - rho[node_cells[0]]) / Vj[node_cells[0]];
if (node_cells.size() == 1) { Gradj[node_cells[0]][1] = (rhoU_node_mean[node_id][0] - rhoU[node_cells[0]][0]) / Vj[node_cells[0]];
continue; Gradj[node_cells[0]][2] = (rhoE_node_mean[node_id] - rhoE[node_cells[0]]) / Vj[node_cells[0]];
}
if (l == 0) {
Gradj[node_cells[l]][0] = (rho_node_mean[node_id] - rho[node_cells[l]]) / Vj[node_cells[l]];
Gradj[node_cells[l]][1] = (rhoU_node_mean[node_id][0] - rhoU[node_cells[l]][0]) / Vj[node_cells[l]];
Gradj[node_cells[l]][2] = (rhoE_node_mean[node_id] - rhoE[node_cells[l]]) / Vj[node_cells[l]];
} else {
Gradjpu[node_cells[l]][0] = (rho[node_cells[l]] - rho_node_mean[node_id]) / Vj[node_cells[l]];
Gradjpu[node_cells[l]][1] = (rhoU[node_cells[l]][0] - rhoU_node_mean[node_id][0]) / Vj[node_cells[l]];
Gradjpu[node_cells[l]][2] = (rhoE[node_cells[l]] - rhoE_node_mean[node_id]) / Vj[node_cells[l]];
;
}
}
}); });
return {Gradj, Gradjpu}; return Gradj;
} }
TinyVector<SIZEproblem> TinyVector<SIZEproblem>
...@@ -247,7 +235,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -247,7 +235,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
return B; return B;
} }
std::tuple<CellValue<TinyVector<SIZEproblem>>, CellValue<TinyVector<SIZEproblem>>> CellValue<TinyVector<SIZEproblem>>
_computeGradConservVariables_for_b(const MeshType& mesh, _computeGradConservVariables_for_b(const MeshType& mesh,
const DiscreteScalarFunction& rho, const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& rhoU, const DiscreteVectorFunction& rhoU,
...@@ -264,31 +252,17 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -264,31 +252,17 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()}; CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()};
Gradj.fill(zero); Gradj.fill(zero);
CellValue<TinyVector<SIZEproblem>> Gradjpu{mesh.connectivity()};
Gradjpu.fill(zero);
parallel_for( parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id]; const auto& node_cells = node_to_cell_matrix[node_id];
for (size_t l = 0; l < node_cells.size(); l++) { Gradj[node_cells[0]][0] = (rho_node[node_id] - rho[node_cells[0]]) / Vj[node_cells[0]];
if (node_cells.size() == 1) { Gradj[node_cells[0]][1] = (rhoU_node[node_id][0] - rhoU[node_cells[0]][0]) / Vj[node_cells[0]];
continue; Gradj[node_cells[0]][2] = (rhoE_node[node_id] - rhoE[node_cells[0]]) / Vj[node_cells[0]];
}
if (l == 0) {
Gradj[node_cells[l]][0] = (rho_node[node_id] - rho[node_cells[l]]) / Vj[node_cells[l]];
Gradj[node_cells[l]][1] = (rhoU_node[node_id][0] - rhoU[node_cells[l]][0]) / Vj[node_cells[l]];
Gradj[node_cells[l]][2] = (rhoE_node[node_id] - rhoE[node_cells[l]]) / Vj[node_cells[l]];
} else {
Gradjpu[node_cells[l]][0] = (rho[node_cells[l]] - rho_node[node_id]) / Vj[node_cells[l]];
Gradjpu[node_cells[l]][1] = (rhoU[node_cells[l]][0] - rhoU_node[node_id][0]) / Vj[node_cells[l]];
Gradjpu[node_cells[l]][2] = (rhoE[node_cells[l]] - rhoE_node[node_id]) / Vj[node_cells[l]];
;
}
}
}); });
return {Gradj, Gradjpu}; return Gradj;
} }
std::tuple<TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>> std::tuple<TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>>
...@@ -444,41 +418,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -444,41 +418,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
double U_2 = U[cell_id][0] * U[cell_id][0];
err_function_term[cell_id] = std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]); err_function_term[cell_id] = 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]); exp_term[cell_id] =
std::exp(-lambda[cell_id] * U[cell_id][0] * U[cell_id][0]) / std::sqrt(pi * lambda[cell_id]);
}); });
// NodeValue<double> rho_node_mean{mesh.connectivity()};
// NodeValue<Rd> rhoU_node_mean{mesh.connectivity()};
// NodeValue<double> rhoE_node_mean{mesh.connectivity()};
// NodeValue<Rd> U_node_mean{mesh.connectivity()};
// NodeValue<double> lambda_node_mean{mesh.connectivity()};
// rho_node_mean.fill(1);
// rhoU_node_mean.fill(zero);
// rhoE_node_mean.fill(1);
// U_node_mean.fill(zero);
// lambda_node_mean.fill(1);
// parallel_for(
// mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
// const auto& node_cells = node_to_cell_matrix[node_id];
// for (size_t l = 0; l < node_cells.size(); l++) {
// if (node_cells.size() == 1)
// continue;
// rho_node_mean[node_id] += rho[node_cells[l]];
// rhoU_node_mean[node_id] += rhoU[node_cells[l]];
// rhoE_node_mean[node_id] += rhoE[node_cells[l]];
// }
// rho_node_mean[node_id] *= 0.5;
// rhoU_node_mean[node_id] *= 0.5;
// rhoE_node_mean[node_id] *= 0.5;
// });
parallel_for( parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id]; const auto& node_cells = node_to_cell_matrix[node_id];
...@@ -541,9 +485,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -541,9 +485,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
} }
} }
auto [Gradj_for_a, Gradjpu_for_a] = this->_computeGradConservVariables_for_a(mesh, rho, rhoU, rhoE); auto Gradj_for_a = this->_computeGradConservVariables_for_a(mesh, rho, rhoU, rhoE);
auto [Gradj_for_b, Gradjpu_for_b] = auto Gradj_for_b = this->_computeGradConservVariables_for_b(mesh, rho, rhoU, rhoE, rho_node, rhoU_node, rhoE_node);
this->_computeGradConservVariables_for_b(mesh, rho, rhoU, rhoE, rho_node, rhoU_node, rhoE_node);
NodeValue<TinyVector<SIZEproblem>> Fluxes_term{mesh.connectivity()}; NodeValue<TinyVector<SIZEproblem>> Fluxes_term{mesh.connectivity()};
Fluxes_term.fill(zero); Fluxes_term.fill(zero);
...@@ -567,8 +510,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -567,8 +510,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoU_right_node[0] = rhoU[node_cells[1]][0] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][1]; rhoU_right_node[0] = rhoU[node_cells[1]][0] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][1];
const double rhoE_right_node = rhoE[node_cells[1]] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][2]; const double rhoE_right_node = rhoE[node_cells[1]] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][2];
// peut etre des + 0.5* ici (au dessus)
Rd U_right_node; Rd U_right_node;
U_right_node[0] = rhoU_right_node[0] / rho_right_node; U_right_node[0] = rhoU_right_node[0] / rho_right_node;
const double lambda_right_node = const double lambda_right_node =
...@@ -580,45 +521,44 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -580,45 +521,44 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
TinyVector<7> u_moments_right = this->_compute_u_moments(rho_right_node, U_right_node, lambda_right_node); TinyVector<7> u_moments_right = this->_compute_u_moments(rho_right_node, U_right_node, lambda_right_node);
TinyVector<2> xi2_moments_right = this->_compute_xi2_moments(lambda_right_node, delta); TinyVector<2> xi2_moments_right = this->_compute_xi2_moments(lambda_right_node, delta);
// TinyMatrix<SIZEproblem> InvM_left =
// this->_computeInvMatrixM(rho_left_node, U_left_node, lambda_left_node, delta);
// TinyMatrix<SIZEproblem> InvM_right =
// this->_computeInvMatrixM(rho_right_node, U_right_node, lambda_right_node, delta);
// TinyMatrix<SIZEproblem> InvMInterface =
// this->_computeInvMatrixM(rho_node[node_id], U_node[node_id], lambda_node[node_id], delta);
TinyMatrix<SIZEproblem> InvM_left = TinyMatrix<SIZEproblem> InvM_left =
this->_computeInvM_bis(rho_left_node, rhoU_left_node, rhoE_left_node, U_left_node, lambda_left_node, delta); this->_computeInvMatrixM(rho_left_node, U_left_node, lambda_left_node, delta);
TinyMatrix<SIZEproblem> InvM_right = this->_computeInvM_bis(rho_right_node, rhoU_right_node, rhoE_right_node, TinyMatrix<SIZEproblem> InvM_right =
U_right_node, lambda_right_node, delta); this->_computeInvMatrixM(rho_right_node, U_right_node, lambda_right_node, delta);
TinyMatrix<SIZEproblem> InvMInterface = TinyMatrix<SIZEproblem> InvMInterface =
this->_computeInvM_bis(rho_node[node_id], rhoU_node[node_id], rhoE_node[node_id], U_node[node_id], this->_computeInvMatrixM(rho_node[node_id], U_node[node_id], lambda_node[node_id], delta);
lambda_node[node_id], delta);
// std::cout << std::endl << "InvM : " << InvM << "\t and \t InvM_bis" << InvM_bis << std::endl; // TinyMatrix<SIZEproblem> InvM_left =
// this->_computeInvM_bis(rho_left_node, rhoU_left_node, rhoE_left_node, U_left_node, lambda_left_node,
// delta);
// TinyMatrix<SIZEproblem> InvM_right = this->_computeInvM_bis(rho_right_node, rhoU_right_node, rhoE_right_node,
// U_right_node, lambda_right_node, delta);
// TinyMatrix<SIZEproblem> InvMInterface =
// this->_computeInvM_bis(rho_node[node_id], rhoU_node[node_id], rhoE_node[node_id], U_node[node_id],
// lambda_node[node_id], delta);
TinyVector<SIZEproblem> al; TinyVector<SIZEproblem> al;
TinyVector<SIZEproblem> ar; TinyVector<SIZEproblem> ar;
TinyVector<SIZEproblem> bl; TinyVector<SIZEproblem> bl;
TinyVector<SIZEproblem> br; TinyVector<SIZEproblem> br;
for (size_t l = 0; l < node_cells.size(); l++) { al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[0]]);
if (node_cells.size() == 1) { bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[0]]);
continue; ar = this->_compute_a_and_b(InvM_right, Gradj_for_a[node_cells[0]]);
} br = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[0]]);
if (l == 0) {
al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[l]]);
bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[l]]);
} else {
ar = this->_compute_a_and_b(InvM_right, Gradjpu_for_a[node_cells[l]]);
br = this->_compute_a_and_b(InvMInterface, Gradjpu_for_b[node_cells[l]]);
}
}
// TinyMatrix<SIZEproblem> Matrix_for_A_left = // for (size_t l = 0; l < node_cells.size(); l++) {
// this->_computeMatricesC1_full_moments(u_moments_left, xi2_moments_left); // if (node_cells.size() == 1) {
// TinyMatrix<SIZEproblem> Matrix_for_A_right = // continue;
// this->_computeMatricesC1_full_moments(u_moments_right, xi2_moments_right); // }
// if (l == 0) {
// al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[l]]);
// bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[l]]);
// } else {
// ar = this->_compute_a_and_b(InvM_right, -Gradj_for_a[node_cells[l]]);
// br = this->_compute_a_and_b(InvMInterface, -Gradj_for_b[node_cells[l]]);
// }
// }
auto [C0_semi_pos_left, C1_semi_pos_left, C2_semi_pos_left] = auto [C0_semi_pos_left, C1_semi_pos_left, C2_semi_pos_left] =
this->_computeMatricesC0C1C2_semi_moments(rho_left_node, U_left_node, lambda_left_node, u_moments_left, this->_computeMatricesC0C1C2_semi_moments(rho_left_node, U_left_node, lambda_left_node, u_moments_left,
...@@ -668,8 +608,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER ...@@ -668,8 +608,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
F_fn[1] = C1_semi_pos_left(1, 0) + C1_semi_neg_right(1, 0); F_fn[1] = C1_semi_pos_left(1, 0) + C1_semi_neg_right(1, 0);
F_fn[2] = C1_semi_pos_left(2, 0) + C1_semi_neg_right(2, 0); F_fn[2] = C1_semi_pos_left(2, 0) + C1_semi_neg_right(2, 0);
Fluxes_term[node_id] = (1 - eta[node_id]) * G + eta[node_id] * F_fn + Fluxes_from_a + Fluxes_from_b + Fluxes_term[node_id] = (1 - eta[node_id]) * G + eta[node_id] * F_fn + Fluxes_for_Navier + Fluxes_from_a +
Fluxes_from_B + Fluxes_for_Navier; Fluxes_from_b + Fluxes_from_B;
}); });
NodeValuePerCell<TinyVector<SIZEproblem>> Fluxes_term_apply{mesh.connectivity()}; NodeValuePerCell<TinyVector<SIZEproblem>> Fluxes_term_apply{mesh.connectivity()};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment