petsc_diff_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
23 #include "libmesh/petsc_matrix.h"
24 #include "libmesh/petsc_vector.h"
26 
27 #ifdef LIBMESH_HAVE_PETSC
28 
29 namespace libMesh
30 {
31 
32 //--------------------------------------------------------------------
33 // Functions with C linkage to pass to PETSc. PETSc will call these
34 // methods as needed.
35 //
36 // Since they must have C linkage they have no knowledge of a namespace.
37 // Give them an obscure name to avoid namespace pollution.
38 extern "C"
39 {
40  // Function to hand to PETSc's SNES,
41  // which monitors convergence at X
42  PetscErrorCode
44  PetscInt its,
45  PetscReal fnorm,
46  void * ctx)
47  {
48  PetscDiffSolver & solver =
49  *(static_cast<PetscDiffSolver *> (ctx));
50 
51  if (solver.verbose)
52  libMesh::out << " PetscDiffSolver step " << its
53  << ", |residual|_2 = " << fnorm << std::endl;
54  if (solver.linear_solution_monitor.get()) {
55  int ierr = 0;
56 
57  Vec petsc_delta_u;
58  ierr = SNESGetSolutionUpdate(snes, &petsc_delta_u);
59  CHKERRABORT(solver.comm().get(), ierr);
60  PetscVector<Number> delta_u(petsc_delta_u, solver.comm());
61  delta_u.close();
62 
63  Vec petsc_u;
64  ierr = SNESGetSolution(snes, &petsc_u);
65  CHKERRABORT(solver.comm().get(), ierr);
66  PetscVector<Number> u(petsc_u, solver.comm());
67  u.close();
68 
69  Vec petsc_res;
70  ierr = SNESGetFunction(snes, &petsc_res, nullptr, nullptr);
71  CHKERRABORT(solver.comm().get(), ierr);
72  PetscVector<Number> res(petsc_res, solver.comm());
73  res.close();
74 
75  (*solver.linear_solution_monitor)(
76  delta_u, delta_u.l2_norm(),
77  u, u.l2_norm(),
78  res, res.l2_norm(), its);
79  }
80  return 0;
81  }
82 
83  // Functions to hand to PETSc's SNES,
84  // which compute the residual or jacobian at X
85  PetscErrorCode
86  __libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void * ctx)
87  {
88  libmesh_assert(x);
89  libmesh_assert(r);
90  libmesh_assert(ctx);
91 
92  PetscDiffSolver & solver =
93  *(static_cast<PetscDiffSolver*> (ctx));
94  ImplicitSystem & sys = solver.system();
95 
96  if (solver.verbose)
97  libMesh::out << "Assembling the residual" << std::endl;
98 
99  PetscVector<Number> & X_system =
100  *cast_ptr<PetscVector<Number> *>(sys.solution.get());
101  PetscVector<Number> & R_system =
102  *cast_ptr<PetscVector<Number> *>(sys.rhs);
103  PetscVector<Number> X_input(x, sys.comm()), R_input(r, sys.comm());
104 
105  // DiffSystem assembles from the solution and into the rhs, so swap
106  // those with our input vectors before assembling. They'll probably
107  // already be references to the same vectors, but PETSc might do
108  // something tricky.
109  X_input.swap(X_system);
110  R_input.swap(R_system);
111 
112  // We may need to localize a parallel solution
113  sys.update();
114 
115  // We may need to correct a non-conforming solution
117 
118  // Do DiffSystem assembly
119  sys.assembly(true, false);
120  R_system.close();
121 
122  // Swap back
123  X_input.swap(X_system);
124  R_input.swap(R_system);
125 
126  // No errors, we hope
127  return 0;
128  }
129 
130 
131  PetscErrorCode
133  Vec x,
134 #if PETSC_RELEASE_LESS_THAN(3,5,0)
135  Mat * libmesh_dbg_var(j),
136  Mat * pc,
137  MatStructure * msflag,
138 #else
139  Mat libmesh_dbg_var(j),
140  Mat pc,
141 #endif
142  void * ctx)
143  {
144  libmesh_assert(x);
145  libmesh_assert(j);
146  // libmesh_assert_equal_to (pc, j); // We don't use separate preconditioners yet
147  libmesh_assert(ctx);
148 
149  PetscDiffSolver & solver =
150  *(static_cast<PetscDiffSolver*> (ctx));
151  ImplicitSystem & sys = solver.system();
152 
153  if (solver.verbose)
154  libMesh::out << "Assembling the Jacobian" << std::endl;
155 
156  PetscVector<Number> & X_system =
157  *cast_ptr<PetscVector<Number> *>(sys.solution.get());
158  PetscVector<Number> X_input(x, sys.comm());
159 
160 #if PETSC_RELEASE_LESS_THAN(3,5,0)
161  PetscMatrix<Number> J_input(*pc, sys.comm());
162 #else
163  PetscMatrix<Number> J_input(pc, sys.comm());
164 #endif
165  PetscMatrix<Number> & J_system =
166  *cast_ptr<PetscMatrix<Number> *>(sys.matrix);
167 
168  // DiffSystem assembles from the solution and into the jacobian, so
169  // swap those with our input vectors before assembling. They'll
170  // probably already be references to the same vectors, but PETSc
171  // might do something tricky.
172  X_input.swap(X_system);
173  J_input.swap(J_system);
174 
175  // We may need to localize a parallel solution
176  sys.update();
177 
178  // We may need to correct a non-conforming solution
180 
181  // Do DiffSystem assembly
182  sys.assembly(false, true);
183  J_system.close();
184 
185  // Swap back
186  X_input.swap(X_system);
187  J_input.swap(J_system);
188 
189 #if PETSC_RELEASE_LESS_THAN(3,5,0)
190  *msflag = SAME_NONZERO_PATTERN;
191 #endif
192  // No errors, we hope
193  return 0;
194  }
195 
196 } // extern "C"
197 
198 
200  : Parent(s)
201 {
202 }
203 
204 
206 {
207  LOG_SCOPE("init()", "PetscDiffSolver");
208 
209  Parent::init();
210 
211  this->setup_petsc_data();
212 }
213 
214 
215 
217 {
218  this->clear();
219 }
220 
221 
222 
224 {
225  LOG_SCOPE("clear()", "PetscDiffSolver");
226 
227  int ierr = LibMeshSNESDestroy(&_snes);
228  LIBMESH_CHKERR(ierr);
229 
230  _dm_wrapper.clear();
231 }
232 
233 
234 
236 {
237  LOG_SCOPE("reinit()", "PetscDiffSolver");
238 
239  // We need to wipe out all the old PETSc data
240  // if we are reinit'ing, since we'll need to build
241  // it all back up again.
242  this->clear();
243 
244  Parent::reinit();
245 
246  this->setup_petsc_data();
247 }
248 
249 
250 
252 {
253  switch (r)
254  {
255  case SNES_CONVERGED_FNORM_ABS:
257  case SNES_CONVERGED_FNORM_RELATIVE:
259 #if PETSC_VERSION_LESS_THAN(3,2,1)
260  case SNES_CONVERGED_PNORM_RELATIVE:
261 #else
262  case SNES_CONVERGED_SNORM_RELATIVE:
263 #endif
265  case SNES_CONVERGED_ITS:
266  case SNES_CONVERGED_TR_DELTA:
268  case SNES_DIVERGED_FUNCTION_DOMAIN:
269  case SNES_DIVERGED_FUNCTION_COUNT:
270  case SNES_DIVERGED_FNORM_NAN:
271 #if !PETSC_VERSION_LESS_THAN(3,3,0)
272  case SNES_DIVERGED_INNER:
273 #endif
274  case SNES_DIVERGED_LINEAR_SOLVE:
275  case SNES_DIVERGED_LOCAL_MIN:
277  case SNES_DIVERGED_MAX_IT:
279 #if PETSC_VERSION_LESS_THAN(3,2,0)
280  case SNES_DIVERGED_LS_FAILURE:
281 #else
282  case SNES_DIVERGED_LINE_SEARCH:
283 #endif
285  // In PETSc, SNES_CONVERGED_ITERATING means
286  // the solve is still iterating, but by the
287  // time we get here, we must have either
288  // converged or diverged, so
289  // SNES_CONVERGED_ITERATING is invalid.
290  case SNES_CONVERGED_ITERATING:
292  default:
293  break;
294  }
296 }
297 
298 
299 
301 {
302  LOG_SCOPE("solve()", "PetscDiffSolver");
303 
304  PetscVector<Number> & x =
305  *(cast_ptr<PetscVector<Number> *>(_system.solution.get()));
306  PetscMatrix<Number> & jac =
307  *(cast_ptr<PetscMatrix<Number> *>(_system.matrix));
308  PetscVector<Number> & r =
309  *(cast_ptr<PetscVector<Number> *>(_system.rhs));
310 
311  int ierr = 0;
312 
313  ierr = SNESSetFunction (_snes, r.vec(),
315  LIBMESH_CHKERR(ierr);
316 
317  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
319  LIBMESH_CHKERR(ierr);
320 
321  ierr = SNESSetFromOptions(_snes);
322  LIBMESH_CHKERR(ierr);
323 
324  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
325  LIBMESH_CHKERR(ierr);
326 
327 #ifdef LIBMESH_ENABLE_CONSTRAINTS
329 #endif
330 
331  SNESConvergedReason reason;
332  SNESGetConvergedReason(_snes, &reason);
333 
334  return convert_solve_result(reason);
335 }
336 
338 {
339  int ierr=0;
340 
341  ierr = SNESCreate(this->comm().get(),&_snes);
342  LIBMESH_CHKERR(ierr);
343 
345  this, PETSC_NULL);
346  LIBMESH_CHKERR(ierr);
347 
348  if (libMesh::on_command_line("--solver-system-names"))
349  {
350  ierr = SNESSetOptionsPrefix(_snes, (_system.name()+"_").c_str());
351  LIBMESH_CHKERR(ierr);
352  }
353 
354  bool use_petsc_dm = libMesh::on_command_line("--use_petsc_dm");
355 
356  // This needs to be called before SNESSetFromOptions
357  if (use_petsc_dm)
359 
360  // If we're not using PETSc DM, let's keep around
361  // the old style for fieldsplit
362  if (!use_petsc_dm)
363  {
364  ierr = SNESSetFromOptions(_snes);
365  LIBMESH_CHKERR(ierr);
366 
367  KSP my_ksp;
368  ierr = SNESGetKSP(_snes, &my_ksp);
369  LIBMESH_CHKERR(ierr);
370 
371  PC my_pc;
372  ierr = KSPGetPC(my_ksp, &my_pc);
373  LIBMESH_CHKERR(ierr);
374 
376  }
377 }
378 
379 } // namespace libMesh
380 
381 #endif // LIBMESH_HAVE_PETSC
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
virtual void init() override
DiffSolver::SolveResult convert_solve_result(SNESConvergedReason r)
NumericVector< Number > * rhs
const Parallel::Communicator & comm() const
virtual void reinit()
Definition: diff_solver.C:60
virtual unsigned int solve() override
void swap(PetscMatrix< T > &)
PetscDiffSolver(sys_type &system)
virtual void swap(NumericVector< T > &v) override
virtual void init()
Definition: diff_solver.C:69
PetscErrorCode __libmesh_petsc_diff_solver_jacobian(SNES, Vec x, #if PETSC_RELEASE_LESS_THAN(3, 5, 0) Mat *libmesh_dbg_var(j), Mat *pc, MatStructure *msflag, #else Mat libmesh_dbg_var(j), Mat pc, #endif void *ctx)
virtual void assembly(bool, bool, bool=false, bool=false)
sys_type & _system
Definition: diff_solver.h:319
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
const sys_type & system() const
Definition: diff_solver.h:138
virtual void update()
Definition: system.C:408
PetscErrorCode ierr
SparseMatrix< Number > * matrix
SparseMatrix interface to PETSc Mat.
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1535
bool on_command_line(std::string arg)
Definition: libmesh.C:876
const std::string & name() const
Definition: system.h:2017
PetscErrorCode __libmesh_petsc_diff_solver_residual(SNES, Vec x, Vec r, void *ctx)
OStreamProxy out(std::cout)
void clear()
Destroys and clears all build DM-related data.
const DofMap & get_dof_map() const
Definition: system.h:2049
virtual void reinit() override
void init_and_attach_petscdm(System &system, SNES &snes)
virtual void close() override
Definition: petsc_vector.h:810
std::unique_ptr< LinearSolutionMonitor > linear_solution_monitor
Definition: diff_solver.h:289
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
PetscErrorCode __libmesh_petsc_diff_solver_monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
Manages consistently variables, degrees of freedom, coefficient vectors, and matrices for implicit sy...