From 1cf9fa3b7e213e167e8f5dcc29b39aa774ba0d23 Mon Sep 17 00:00:00 2001
From: maxime <roussel.maxime012@gmail.com>
Date: Fri, 26 Jul 2024 09:29:29 +0200
Subject: [PATCH] premier schema
---
src/scheme/ImplicitAcousticO2Solver.cpp | 41 +++++++++++++++++++++----
1 file changed, 35 insertions(+), 6 deletions(-)
diff --git a/src/scheme/ImplicitAcousticO2Solver.cpp b/src/scheme/ImplicitAcousticO2Solver.cpp
index 9351d9119..42af7235f 100644
--- a/src/scheme/ImplicitAcousticO2Solver.cpp
+++ b/src/scheme/ImplicitAcousticO2Solver.cpp
@@ -1707,15 +1707,44 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
int number_iter = 0;
do {
number_iter++;
+ // first step
- this->_computeGradJU_AU(u, p, pi, dt);
+ this->_computeGradJU_AU(u, p, pi, 0.5 * dt);
NodeValue<const Rd> ur = this->_getUr();
NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
+ // final step
- // time n+1
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ CellValue<double> final_step_tau(m_mesh.connectivity());
+ CellValue<Rd> final_step_u = copy(u.cellValues());
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ Rd momentum_fluxes = zero;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ momentum_fluxes += Fjr(j, R);
+ }
+ const double dt_over_Mj = dt / m_Mj[j];
+ final_step_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
+ });
+
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ double final_step_tau_fluxes = 0;
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ final_step_tau_fluxes += dot(m_Djr(j, R), ur[r]);
+ }
+ final_step_tau[j] = m_tau[j] + 0.5 * (dt / m_Mj[j]) * final_step_tau_fluxes;
+ };
+
+ // time n+1
+
NodeValue<Rd> new_xr = copy(mesh->xr());
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
@@ -1725,7 +1754,7 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
new_rho = copy(rho.cellValues());
- new_u = copy(u.cellValues());
+ new_u = final_step_u;
new_E = copy(E.cellValues());
parallel_for(
@@ -1740,8 +1769,8 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const double dt_over_Mj = dt / m_Mj[j];
- new_u[j] -= dt_over_Mj * momentum_fluxes;
- new_E[j] -= dt_over_Mj * energy_fluxes;
+ new_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
+ new_E[j] -= 0.5 * dt_over_Mj * energy_fluxes;
});
// update Djr
@@ -1786,7 +1815,7 @@ class ImplicitAcousticO2SolverHandler::ImplicitAcousticO2Solver final
const NodeId r = cell_nodes[R];
new_tau_fluxes += dot(m_Djr(j, R), ur[r]);
}
- new_tau[j] = m_tau[j] + (dt / m_Mj[j]) * new_tau_fluxes;
+ new_tau[j] = final_step_tau[j] + 0.5 * (dt / m_Mj[j]) * new_tau_fluxes;
// std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
}
--
GitLab