Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/diff_system.h"
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/libmesh_logging.h"
22 : #include "libmesh/petsc_diff_solver.h"
23 : #include "libmesh/petsc_matrix_base.h"
24 : #include "libmesh/petsc_vector.h"
25 : #include "libmesh/petsc_auto_fieldsplit.h"
26 : #include "libmesh/boundary_info.h"
27 :
28 : #ifdef LIBMESH_HAVE_PETSC
29 :
30 : namespace libMesh
31 : {
32 :
33 : //--------------------------------------------------------------------
34 : // Functions with C linkage to pass to PETSc. PETSc will call these
35 : // methods as needed.
36 : //
37 : // Since they must have C linkage they have no knowledge of a namespace.
38 : // Give them an obscure name to avoid namespace pollution.
39 : extern "C"
40 : {
41 : // Function to hand to PETSc's SNES,
42 : // which monitors convergence at X
43 : PetscErrorCode
44 9890 : __libmesh_petsc_diff_solver_monitor (SNES snes,
45 : PetscInt its,
46 : PetscReal fnorm,
47 : void * ctx)
48 : {
49 : PetscFunctionBegin;
50 :
51 280 : PetscDiffSolver & solver =
52 : *(static_cast<PetscDiffSolver *> (ctx));
53 :
54 9890 : if (solver.verbose)
55 120 : libMesh::out << " PetscDiffSolver step " << its
56 120 : << ", |residual|_2 = " << fnorm << std::endl;
57 9890 : if (solver.linear_solution_monitor.get())
58 : {
59 : Vec petsc_delta_u;
60 0 : LibmeshPetscCall2(solver.comm(), SNESGetSolutionUpdate(snes, &petsc_delta_u));
61 0 : PetscVector<Number> delta_u(petsc_delta_u, solver.comm());
62 0 : delta_u.close();
63 :
64 : Vec petsc_u;
65 0 : LibmeshPetscCall2(solver.comm(), SNESGetSolution(snes, &petsc_u));
66 0 : PetscVector<Number> u(petsc_u, solver.comm());
67 0 : u.close();
68 :
69 : Vec petsc_res;
70 0 : LibmeshPetscCall2(solver.comm(), SNESGetFunction(snes, &petsc_res, nullptr, nullptr));
71 0 : PetscVector<Number> res(petsc_res, solver.comm());
72 0 : res.close();
73 :
74 0 : (*solver.linear_solution_monitor)(
75 0 : delta_u, delta_u.l2_norm(),
76 0 : u, u.l2_norm(),
77 0 : res, res.l2_norm(), its);
78 0 : }
79 9890 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
80 : }
81 :
82 : // Functions to hand to PETSc's SNES,
83 : // which compute the residual or jacobian at X
84 : PetscErrorCode
85 12455 : __libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void * ctx)
86 : {
87 : PetscFunctionBegin;
88 :
89 352 : libmesh_assert(x);
90 352 : libmesh_assert(r);
91 352 : libmesh_assert(ctx);
92 :
93 352 : PetscDiffSolver & solver =
94 : *(static_cast<PetscDiffSolver*> (ctx));
95 704 : ImplicitSystem & sys = solver.system();
96 :
97 12455 : if (solver.verbose)
98 192 : libMesh::out << "Assembling the residual" << std::endl;
99 :
100 : PetscVector<Number> & X_system =
101 352 : *cast_ptr<PetscVector<Number> *>(sys.solution.get());
102 : PetscVector<Number> & R_system =
103 12455 : *cast_ptr<PetscVector<Number> *>(sys.rhs);
104 13159 : PetscVector<Number> X_input(x, sys.comm()), R_input(r, sys.comm());
105 :
106 : // DiffSystem assembles from the solution and into the rhs, so swap
107 : // those with our input vectors before assembling. They'll probably
108 : // already be references to the same vectors, but PETSc might do
109 : // something tricky.
110 12455 : X_input.swap(X_system);
111 12455 : R_input.swap(R_system);
112 :
113 : // We may need to localize a parallel solution
114 12455 : sys.update();
115 :
116 : // We may need to correct a non-conforming solution
117 12455 : if (solver.exact_constraint_enforcement())
118 12455 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
119 :
120 : // Do DiffSystem assembly
121 12455 : sys.assembly(true, false, !solver.exact_constraint_enforcement());
122 12455 : R_system.close();
123 :
124 : // Swap back
125 12455 : X_input.swap(X_system);
126 12455 : R_input.swap(R_system);
127 :
128 : // No errors, we hope
129 12807 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
130 11751 : }
131 :
132 :
133 : PetscErrorCode
134 6670 : __libmesh_petsc_diff_solver_jacobian (SNES,
135 : Vec x,
136 : Mat libmesh_dbg_var(j),
137 : Mat libmesh_dbg_var(pc),
138 : void * ctx)
139 : {
140 : PetscFunctionBegin;
141 :
142 188 : libmesh_assert(x);
143 188 : libmesh_assert(j);
144 : // libmesh_assert_equal_to (pc, j); // We don't use separate preconditioners yet
145 188 : libmesh_assert(ctx);
146 :
147 188 : PetscDiffSolver & solver =
148 : *(static_cast<PetscDiffSolver*> (ctx));
149 376 : ImplicitSystem & sys = solver.system();
150 :
151 6670 : if (solver.verbose)
152 108 : libMesh::out << "Assembling the Jacobian" << std::endl;
153 :
154 : PetscVector<Number> & X_system =
155 188 : *cast_ptr<PetscVector<Number> *>(sys.solution.get());
156 6858 : PetscVector<Number> X_input(x, sys.comm());
157 :
158 : PetscMatrixBase<Number> & J_system =
159 6670 : *cast_ptr<PetscMatrixBase<Number> *>(sys.matrix);
160 188 : libmesh_assert(J_system.mat() == pc);
161 :
162 : // DiffSystem assembles from the solution and into the jacobian, so
163 : // swap those with our input vectors before assembling. They'll
164 : // probably already be references to the same vectors, but PETSc
165 : // might do something tricky.
166 6670 : X_input.swap(X_system);
167 :
168 : // We may need to localize a parallel solution
169 6670 : sys.update();
170 :
171 : // We may need to correct a non-conforming solution
172 6670 : if (solver.exact_constraint_enforcement())
173 6670 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
174 :
175 : // Do DiffSystem assembly
176 6670 : sys.assembly(false, true, !solver.exact_constraint_enforcement());
177 6670 : J_system.close();
178 :
179 : // Swap back
180 6670 : X_input.swap(X_system);
181 :
182 : // No errors, we hope
183 6858 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
184 6294 : }
185 :
186 : } // extern "C"
187 :
188 :
189 420 : PetscDiffSolver::PetscDiffSolver (sys_type & s)
190 420 : : Parent(s)
191 : {
192 420 : }
193 :
194 :
195 420 : void PetscDiffSolver::init ()
196 : {
197 24 : LOG_SCOPE("init()", "PetscDiffSolver");
198 :
199 420 : Parent::init();
200 :
201 420 : this->setup_petsc_data();
202 420 : }
203 :
204 :
205 :
206 792 : PetscDiffSolver::~PetscDiffSolver () = default;
207 :
208 :
209 :
210 0 : void PetscDiffSolver::clear()
211 : {
212 0 : LOG_SCOPE("clear()", "PetscDiffSolver");
213 :
214 : // calls custom deleter
215 0 : _snes.destroy();
216 :
217 : #if !PETSC_VERSION_LESS_THAN(3,7,3)
218 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
219 0 : _dm_wrapper.clear();
220 : #endif
221 : #endif
222 0 : }
223 :
224 :
225 :
226 0 : void PetscDiffSolver::reinit()
227 : {
228 0 : LOG_SCOPE("reinit()", "PetscDiffSolver");
229 :
230 : // We need to wipe out all the old PETSc data
231 : // if we are reinit'ing, since we'll need to build
232 : // it all back up again.
233 0 : this->clear();
234 :
235 0 : Parent::reinit();
236 :
237 0 : this->setup_petsc_data();
238 0 : }
239 :
240 :
241 :
242 3220 : DiffSolver::SolveResult convert_solve_result(SNESConvergedReason r)
243 : {
244 3220 : switch (r)
245 : {
246 0 : case SNES_CONVERGED_FNORM_ABS:
247 0 : return DiffSolver::CONVERGED_ABSOLUTE_RESIDUAL;
248 92 : case SNES_CONVERGED_FNORM_RELATIVE:
249 92 : return DiffSolver::CONVERGED_RELATIVE_RESIDUAL;
250 0 : case SNES_CONVERGED_SNORM_RELATIVE:
251 0 : return DiffSolver::CONVERGED_RELATIVE_STEP;
252 0 : case SNES_CONVERGED_ITS:
253 : // SNES_CONVERGED_TR_DELTA was changed to a diverged condition,
254 : // SNES_DIVERGED_TR_DELTA, in PETSc 1c6b2ff8df. This change will
255 : // likely be in 3.12 and later releases.
256 : #if PETSC_VERSION_LESS_THAN(3,12,0)
257 : case SNES_CONVERGED_TR_DELTA:
258 : #endif
259 0 : return DiffSolver::CONVERGED_NO_REASON;
260 0 : case SNES_DIVERGED_FUNCTION_DOMAIN:
261 : case SNES_DIVERGED_FUNCTION_COUNT:
262 : case SNES_DIVERGED_FNORM_NAN:
263 : case SNES_DIVERGED_INNER:
264 : case SNES_DIVERGED_LINEAR_SOLVE:
265 : case SNES_DIVERGED_LOCAL_MIN:
266 0 : return DiffSolver::DIVERGED_NO_REASON;
267 0 : case SNES_DIVERGED_MAX_IT:
268 0 : return DiffSolver::DIVERGED_MAX_NONLINEAR_ITERATIONS;
269 0 : case SNES_DIVERGED_LINE_SEARCH:
270 0 : return DiffSolver::DIVERGED_BACKTRACKING_FAILURE;
271 : // In PETSc, SNES_CONVERGED_ITERATING means
272 : // the solve is still iterating, but by the
273 : // time we get here, we must have either
274 : // converged or diverged, so
275 : // SNES_CONVERGED_ITERATING is invalid.
276 0 : case SNES_CONVERGED_ITERATING:
277 0 : return DiffSolver::INVALID_SOLVE_RESULT;
278 0 : default:
279 0 : break;
280 : }
281 0 : return DiffSolver::INVALID_SOLVE_RESULT;
282 : }
283 :
284 :
285 :
286 3220 : unsigned int PetscDiffSolver::solve()
287 : {
288 184 : LOG_SCOPE("solve()", "PetscDiffSolver");
289 :
290 : #if !PETSC_VERSION_LESS_THAN(3,7,3)
291 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
292 : // GMG is currently not supported if we enable children to be associated with
293 : // boundary sides
294 92 : libmesh_assert(!_system.get_mesh().get_boundary_info().is_children_on_boundary_side());
295 : #endif
296 : #endif
297 :
298 : PetscVector<Number> & x =
299 3220 : *(cast_ptr<PetscVector<Number> *>(_system.solution.get()));
300 : PetscMatrixBase<Number> & jac =
301 3220 : *(cast_ptr<PetscMatrixBase<Number> *>(_system.matrix));
302 : PetscVector<Number> & r =
303 3220 : *(cast_ptr<PetscVector<Number> *>(_system.rhs));
304 :
305 3220 : LibmeshPetscCall(SNESSetFunction (_snes, r.vec(),
306 : __libmesh_petsc_diff_solver_residual, this));
307 :
308 3220 : LibmeshPetscCall(SNESSetJacobian (_snes, jac.mat(), jac.mat(),
309 : __libmesh_petsc_diff_solver_jacobian, this));
310 :
311 3220 : LibmeshPetscCall(SNESSetFromOptions(_snes));
312 :
313 3220 : LibmeshPetscCall(SNESSolve (_snes, LIBMESH_PETSC_NULLPTR, x.vec()));
314 :
315 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
316 3220 : if (this->_exact_constraint_enforcement)
317 3312 : _system.get_dof_map().enforce_constraints_exactly(_system);
318 : #endif
319 :
320 : SNESConvergedReason reason;
321 3220 : LibmeshPetscCall(SNESGetConvergedReason(_snes, &reason));
322 :
323 : PetscInt l_its, nl_its;
324 3220 : LibmeshPetscCall(SNESGetLinearSolveIterations(_snes, &l_its));
325 3220 : this->_inner_iterations = l_its;
326 :
327 3220 : LibmeshPetscCall(SNESGetIterationNumber(_snes, &nl_its));
328 3220 : this->_outer_iterations = nl_its;
329 :
330 3404 : return convert_solve_result(reason);
331 : }
332 :
333 420 : void PetscDiffSolver::setup_petsc_data()
334 : {
335 420 : LibmeshPetscCall(SNESCreate(this->comm().get(), _snes.get()));
336 :
337 420 : LibmeshPetscCall(SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor,
338 : this, LIBMESH_PETSC_NULLPTR));
339 :
340 828 : if (libMesh::on_command_line("--solver-system-names"))
341 0 : LibmeshPetscCall(SNESSetOptionsPrefix(_snes, (_system.name()+"_").c_str()));
342 :
343 420 : bool use_petsc_dm = libMesh::on_command_line("--use_petsc_dm");
344 :
345 : // This needs to be called before SNESSetFromOptions
346 : #if !PETSC_VERSION_LESS_THAN(3,7,3)
347 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
348 420 : if (use_petsc_dm)
349 280 : this->_dm_wrapper.init_and_attach_petscdm(_system, _snes);
350 : #endif
351 : #endif
352 :
353 : // If we're not using PETSc DM, let's keep around
354 : // the old style for fieldsplit
355 24 : if (!use_petsc_dm)
356 : {
357 140 : LibmeshPetscCall(SNESSetFromOptions(_snes));
358 :
359 : KSP my_ksp;
360 140 : LibmeshPetscCall(SNESGetKSP(_snes, &my_ksp));
361 :
362 : PC my_pc;
363 140 : LibmeshPetscCall(KSPGetPC(my_ksp, &my_pc));
364 :
365 140 : petsc_auto_fieldsplit(my_pc, _system);
366 : }
367 420 : }
368 :
369 : } // namespace libMesh
370 :
371 : #endif // LIBMESH_HAVE_PETSC
|
