From 09489387a6499e15ce93904325a64657d4965c94 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 22 Dec 2023 12:36:46 +0100
Subject: [PATCH] Add forcing of some boundary conditions
The idea is to ensure that axis lines, symmetry planes and Dirichlet
values are preserved exactly.
This is done by changing (enforcing) nodal velocity to satisfy
boundary conditions more precisely
---
src/scheme/ImplicitAcousticSolver.cpp | 58 +++++++++++++++++++--------
1 file changed, 41 insertions(+), 17 deletions(-)
diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index 411c4112c..df37a7591 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -1133,25 +1133,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
- // std::cout << "Hess J = " << Hess_J_crs << '\n';
-
CRSMatrix gradient_f = Hess_J_crs - A;
- // Vector<double> x(gradient_f.numberOfRows());
- // Vector<double> y(gradient_f.numberOfRows());
- // double max_err = 0;
- // for (size_t i = 0; i < gradient_f.numberOfRows(); ++i) {
- // x.fill(0);
- // x[i] = 1;
- // for (size_t j = 0; j < gradient_f.numberOfRows(); ++j) {
- // y.fill(0);
- // y[j] = 1;
-
- // max_err = std::max(max_err, std::abs(dot(x, gradient_f * y) - dot(y, gradient_f * x)));
- // }
- // }
- // std::cout << "MAX NON SYM=" << rang::fgB::cyan << max_err << rang::fg::reset << '\n';
-
HJ_A_t.pause();
return gradient_f;
@@ -1628,6 +1611,47 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
parallel_for(
m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = m_inv_Ar[r] * b[r]; });
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+
+ if constexpr ((Dimension > 1) and (std::is_same_v<AxisBoundaryCondition, T>)) {
+ const Rd& t = bc.direction();
+ const Rdxd txt = tensorProduct(t, t);
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node]; // on fixe le sommet r
+ computed_ur[node_id] = txt * computed_ur[node_id];
+ }
+ } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node]; // on fixe le sommet r
+ computed_ur[node_id] = P * computed_ur[node_id];
+ }
+ } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+ const Array<const NodeId>& node_list = bc.nodeList();
+ const Array<const Rd>& value_list = bc.valueList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ computed_ur[node_id] = value_list[i_node];
+ }
+ }
+ },
+ boundary_condition);
+ }
+
return computed_ur;
}
--
GitLab