OGS
PETScNonlinearSolver.cpp
Go to the documentation of this file.
1
11#ifdef USE_PETSC
12
14
15#include <petscmat.h>
16#include <petscvec.h>
17
18#include "BaseLib/RunTime.h"
19
20namespace
21{
23{
26 std::vector<GlobalVector*>& x;
27 std::vector<GlobalVector*> const& x_prev;
30 int const process_id;
31};
32
33void printConvergenceReason(SNES const& snes_solver,
34 SNESConvergedReason const& reason)
35{
36 const char* strreason;
37 SNESGetConvergedReasonString(snes_solver, &strreason);
38 if (reason < 0)
39 {
40 INFO("PETScSNES diverged reason: {:s}.", strreason);
41 }
42 else
43 {
44 INFO("PETScSNES converged reason: {:s}.", strreason);
45 }
46}
47
48PetscErrorCode updateResidual(SNES /*snes*/, Vec x, Vec petsc_r,
49 void* petsc_context)
50{
51 auto context = static_cast<PetscContext*>(petsc_context);
52
53 DBUG("PETScNonlinearSolver: residual callback called.");
54
55 VecCopy(x, context->x[context->process_id]->getRawVector());
56
57 // Assemble in ogs context.
58 BaseLib::RunTime time_assembly;
59 time_assembly.start();
60 context->system->assemble(context->x, context->x_prev, context->process_id);
61
62 INFO("[time] Assembly took {} s.", time_assembly.elapsed());
63 context->system->getResidual(*context->x[context->process_id],
64 *context->x_prev[context->process_id],
65 *context->r);
66 context->r->finalizeAssembly();
67 context->J->finalizeAssembly();
68
69 context->system->getJacobian(*context->J);
70 context->system->applyKnownSolutionsPETScSNES(
71 *context->J, *context->r, *context->x[context->process_id]);
72
73 VecCopy(context->r->getRawVector(), petsc_r);
74
75 return 0;
76}
77
78PetscErrorCode updateJacobian(SNES /*snes*/, Vec /*x*/, Mat J,
79 Mat /*same as J*/, void* petsc_context)
80{
81 DBUG("PETScNonlinearSolver: jacobian callback called.");
82 // Assume the system is already assembled.
83 // Copy the results into petsc vectors.
84
85 auto context = static_cast<PetscContext*>(petsc_context);
86 MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
87 MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
88 MatCopy(context->J->getRawMatrix(), J, DIFFERENT_NONZERO_PATTERN);
89 return 0;
90}
91} // namespace
92
93namespace NumLib
94{
96 GlobalLinearSolver& /*linear_solver*/, int const maxiter,
97 std::string prefix)
98{
99 SNESCreate(PETSC_COMM_WORLD, &_snes_solver);
100 if (!prefix.empty())
101 {
102 prefix = prefix + "_";
103 SNESSetOptionsPrefix(_snes_solver, prefix.c_str());
104 }
105 // force SNESSolve() to take at least one iteration regardless of the
106 // initial residual norm
107 SNESSetForceIteration(_snes_solver, PETSC_TRUE);
108
109 // Set the maximum iterations.
110 PetscReal atol, rtol, stol;
111 PetscInt maxf;
112 SNESGetTolerances(_snes_solver, &atol, &rtol, &stol, nullptr, &maxf);
113 SNESSetTolerances(_snes_solver, atol, rtol, stol, maxiter, maxf);
114
115 SNESSetFromOptions(_snes_solver);
116#ifndef NDEBUG
117 PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_WORLD);
118#endif // NDEBUG
119}
120
127
133
135 std::vector<GlobalVector*> const& /*x*/,
136 std::vector<GlobalVector*> const& /*x_prev*/, int const /*process_id*/)
137{
139 {
140 return;
141 }
142
143 OGS_FATAL(
144 "Non-equilibrium initial residuum is not implemented for the "
145 "PETScNonlinearSolver.");
146}
147
149 std::vector<GlobalVector*>& x,
150 std::vector<GlobalVector*> const& x_prev,
151 std::function<void(
152 int,
153 std::vector<GlobalVector*> const&)> const& /*postIterationCallback*/,
154 int const process_id)
155{
156 DBUG("PETScNonlinearSolver: solve()");
157 using TimeDiscretizedSystem =
160
161 auto* system = static_cast<TimeDiscretizedSystem*>(_equation_system);
162
163 DBUG("PETScNonlinearSolver: create vectors");
164 // r and J on which the ogs assembly operates.
166 system->getMatrixSpecifications(process_id), _jacobian_id);
168 system->getMatrixSpecifications(process_id), _residual_id);
169
170 // Vectors and matrices used by SNES for solutions. These will be locked
171 // during the SNES' solve call.
173 system->getMatrixSpecifications(process_id), _petsc_jacobian_id);
174 BaseLib::RunTime timer_dirichlet;
175 double time_dirichlet = 0.0;
176
177 timer_dirichlet.start();
178 system->computeKnownSolutions(*x[process_id], process_id);
179 system->applyKnownSolutions(*x[process_id]);
180 time_dirichlet += timer_dirichlet.elapsed();
181 INFO("[time] Applying Dirichlet BCs took {} s.", time_dirichlet);
182
184 system->getMatrixSpecifications(process_id), _petsc_residual_id);
186 system->getMatrixSpecifications(process_id), _petsc_x_id);
187 MathLib::LinAlg::copy(*x[process_id], x_snes); // Initial guess.
188
189 PetscContext petsc_context{_equation_system, x, x_prev, &r, &J, process_id};
190
191 DBUG("PETScNonlinearSolver: set function");
192 SNESSetFunction(_snes_solver, r_snes.getRawVector(), updateResidual,
193 &petsc_context);
194
195 DBUG("PETScNonlinearSolver: set jacobian");
196 // The jacobian and the preconditioner matrices point to the same location.
197 SNESSetJacobian(_snes_solver, J_snes.getRawMatrix(), J_snes.getRawMatrix(),
198 updateJacobian, &petsc_context);
199
200 std::unique_ptr<GlobalVector> xl = nullptr;
201 std::unique_ptr<GlobalVector> xu = nullptr;
202
203 SNESType snes_type;
204 SNESGetType(_snes_solver, &snes_type);
205 if ((std::strcmp(snes_type, SNESVINEWTONRSLS) == 0) ||
206 (std::strcmp(snes_type, SNESVINEWTONSSLS) == 0))
207 {
208 // Set optional constraints via callback.
209 DBUG("PETScNonlinearSolver: set constraints");
211 system->getMatrixSpecifications(process_id));
213 system->getMatrixSpecifications(process_id));
214
215 system->updateConstraints(*xl, *xu, process_id);
218
219 SNESVISetVariableBounds(_snes_solver, xl->getRawVector(),
220 xu->getRawVector());
221 }
222
223 DBUG("PETScNonlinearSolver: call SNESSolve");
224 SNESSolve(_snes_solver, nullptr, x_snes.getRawVector());
225
226 SNESConvergedReason reason;
227 SNESGetConvergedReason(_snes_solver, &reason);
228 printConvergenceReason(_snes_solver, reason);
229
230 PetscInt iterations;
231 SNESGetIterationNumber(_snes_solver, &iterations);
232 INFO("PETScSNES used {} iterations.", iterations);
233
234 // Copy back the solution.
235 MathLib::LinAlg::copy(x_snes, *x[process_id]);
236
242
243 return {reason >= 0, iterations};
244}
245
246} // namespace NumLib
247#endif // USE_PETSC
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the RunTime class.
Count the running time.
Definition RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition RunTime.h:42
void start()
Start the timer.
Definition RunTime.h:32
Global vector based on Eigen vector.
Definition EigenVector.h:25
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
ConvergenceCriterion * _convergence_criterion
NonlinearSolverStatus solve(std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::function< void(int, std::vector< GlobalVector * > const &)> const &postIterationCallback, int const process_id) override
void setEquationSystem(System &eq, ConvergenceCriterion &conv_crit)
Set the nonlinear equation system that will be solved.
PETScNonlinearSolver(GlobalLinearSolver &linear_solver, int maxiter, std::string prefix)
std::size_t _jacobian_id
ID of the Jacobian matrix.
void compensateNonEquilibriumInitialResiduum(bool const value)
std::size_t _residual_id
ID of the residual vector.
void calculateNonEquilibriumInitialResiduum(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id) override
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
PetscErrorCode updateResidual(SNES, Vec x, Vec petsc_r, void *petsc_context)
void printConvergenceReason(SNES const &snes_solver, SNESConvergedReason const &reason)
PetscErrorCode updateJacobian(SNES, Vec, Mat J, Mat, void *petsc_context)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.