From bdb16d78cb0fed1d5e6188b59ca30613f43f6607 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 2 Jul 2024 11:08:17 +0200
Subject: [PATCH] [NL] Copy current Newton applyBCs for PETSc-SNES
The system solved by PETSc-SNES is different one and so far
the methods could be reused. With upcoming changes to the
Newton's version of applying Dirichlet-BCs this is not
going to work anymore.
---
NumLib/ODESolver/NonlinearSystem.h | 7 ++++++
NumLib/ODESolver/PETScNonlinearSolver.cpp | 2 +-
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 22 +++++++++++++++++++
NumLib/ODESolver/TimeDiscretizedODESystem.h | 3 +++
4 files changed, 33 insertions(+), 1 deletion(-)
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index b232e988f18..02d2d279380 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -78,6 +78,13 @@ public:
GlobalMatrix& Jac, GlobalVector& res,
GlobalVector& minus_delta_x) const = 0;
+ //! Apply known solutions to the linearized equation system
+ //! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
+ //! \pre computeKnownSolutions() must have been called before.
+ virtual void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac,
+ GlobalVector& res,
+ GlobalVector& x) const = 0;
+
virtual void updateConstraints(GlobalVector& /*lower*/,
GlobalVector& /*upper*/,
int /*process_id*/) = 0;
diff --git a/NumLib/ODESolver/PETScNonlinearSolver.cpp b/NumLib/ODESolver/PETScNonlinearSolver.cpp
index 7eb3475b026..193b0dc6536 100644
--- a/NumLib/ODESolver/PETScNonlinearSolver.cpp
+++ b/NumLib/ODESolver/PETScNonlinearSolver.cpp
@@ -67,7 +67,7 @@ PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
context->J->finalizeAssembly();
context->system->getJacobian(*context->J);
- context->system->applyKnownSolutionsNewton(
+ context->system->applyKnownSolutionsPETScSNES(
*context->J, *context->r, *context->x[context->process_id]);
VecCopy(context->r->getRawVector(), petsc_r);
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index d7e531afdc3..e618773be5c 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -152,6 +152,28 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
}
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>::
+ applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+ GlobalVector& x) const
+{
+ if (!_known_solutions)
+ {
+ return;
+ }
+
+ using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+ std::vector<IndexType> ids;
+ for (auto const& bc : *_known_solutions)
+ {
+ ids.insert(end(ids), begin(bc.ids), end(bc.ids));
+ }
+
+ // For the Newton method the values must be zero
+ std::vector<double> values(ids.size(), 0);
+ MathLib::applyKnownSolution(Jac, res, x, ids, values);
+}
+
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>::
TimeDiscretizedODESystem(const int process_id, ODE& ode,
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index b53f6ca426b..46dd959689f 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -108,6 +108,9 @@ public:
void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
GlobalVector& minus_delta_x) const override;
+ void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
+ GlobalVector& x) const override;
+
void updateConstraints(GlobalVector& lower, GlobalVector& upper,
int const process_id) override
{
--
GitLab