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, libmesh_nullptr, libmesh_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 
231 
232 
234 {
235  LOG_SCOPE("reinit()", "PetscDiffSolver");
236 
237  // We need to wipe out all the old PETSc data
238  // if we are reinit'ing, since we'll need to build
239  // it all back up again.
240  this->clear();
241 
242  Parent::reinit();
243 
244  this->setup_petsc_data();
245 }
246 
247 
248 
250 {
251  switch (r)
252  {
253  case SNES_CONVERGED_FNORM_ABS:
255  case SNES_CONVERGED_FNORM_RELATIVE:
257 #if PETSC_VERSION_LESS_THAN(3,2,1)
258  case SNES_CONVERGED_PNORM_RELATIVE:
259 #else
260  case SNES_CONVERGED_SNORM_RELATIVE:
261 #endif
263  case SNES_CONVERGED_ITS:
264  case SNES_CONVERGED_TR_DELTA:
266  case SNES_DIVERGED_FUNCTION_DOMAIN:
267  case SNES_DIVERGED_FUNCTION_COUNT:
268  case SNES_DIVERGED_FNORM_NAN:
269 #if !PETSC_VERSION_LESS_THAN(3,3,0)
270  case SNES_DIVERGED_INNER:
271 #endif
272  case SNES_DIVERGED_LINEAR_SOLVE:
273  case SNES_DIVERGED_LOCAL_MIN:
275  case SNES_DIVERGED_MAX_IT:
277 #if PETSC_VERSION_LESS_THAN(3,2,0)
278  case SNES_DIVERGED_LS_FAILURE:
279 #else
280  case SNES_DIVERGED_LINE_SEARCH:
281 #endif
283  // In PETSc, SNES_CONVERGED_ITERATING means
284  // the solve is still iterating, but by the
285  // time we get here, we must have either
286  // converged or diverged, so
287  // SNES_CONVERGED_ITERATING is invalid.
288  case SNES_CONVERGED_ITERATING:
290  default:
291  break;
292  }
294 }
295 
296 
297 
299 {
300  LOG_SCOPE("solve()", "PetscDiffSolver");
301 
302  PetscVector<Number> & x =
303  *(cast_ptr<PetscVector<Number> *>(_system.solution.get()));
304  PetscMatrix<Number> & jac =
305  *(cast_ptr<PetscMatrix<Number> *>(_system.matrix));
306  PetscVector<Number> & r =
307  *(cast_ptr<PetscVector<Number> *>(_system.rhs));
308 
309  int ierr = 0;
310 
311  ierr = SNESSetFunction (_snes, r.vec(),
313  LIBMESH_CHKERR(ierr);
314 
315  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
317  LIBMESH_CHKERR(ierr);
318 
319  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
320  LIBMESH_CHKERR(ierr);
321 
322 #ifdef LIBMESH_ENABLE_CONSTRAINTS
324 #endif
325 
326  SNESConvergedReason reason;
327  SNESGetConvergedReason(_snes, &reason);
328 
329  return convert_solve_result(reason);
330 }
331 
333 {
334  int ierr=0;
335 
336  ierr = SNESCreate(this->comm().get(),&_snes);
337  LIBMESH_CHKERR(ierr);
338 
339  ierr = SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor,
340  this, PETSC_NULL);
341  LIBMESH_CHKERR(ierr);
342 
343  if (libMesh::on_command_line("--solver_system_names"))
344  {
345  ierr = SNESSetOptionsPrefix(_snes, (_system.name()+"_").c_str());
346  LIBMESH_CHKERR(ierr);
347  }
348 
349  ierr = SNESSetFromOptions(_snes);
350  LIBMESH_CHKERR(ierr);
351 
352  KSP my_ksp;
353  ierr = SNESGetKSP(_snes, &my_ksp);
354  LIBMESH_CHKERR(ierr);
355 
356  PC my_pc;
357  ierr = KSPGetPC(my_ksp, &my_pc);
358  LIBMESH_CHKERR(ierr);
359 
361 }
362 
363 } // namespace libMesh
364 
365 #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 reinit() libmesh_override
virtual void init() libmesh_override
const class libmesh_nullptr_t libmesh_nullptr
DiffSolver::SolveResult convert_solve_result(SNESConvergedReason r)
NumericVector< Number > * rhs
virtual void reinit()
Definition: diff_solver.C:60
void swap(PetscMatrix< T > &)
const std::string & name() const
Definition: system.h:1998
PetscDiffSolver(sys_type &system)
virtual void init()
Definition: diff_solver.C:69
virtual void assembly(bool, bool, bool=false, bool=false)
const DofMap & get_dof_map() const
Definition: system.h:2030
sys_type & _system
Definition: diff_solver.h:319
const sys_type & system() const
Definition: diff_solver.h:138
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
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)
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
virtual void close() libmesh_override
Definition: petsc_vector.h:807
virtual void update()
Definition: system.C:406
virtual unsigned int solve() libmesh_override
PetscErrorCode ierr
bool on_command_line(const std::string &arg)
Definition: libmesh.C:920
const Parallel::Communicator & comm() const
SparseMatrix< Number > * matrix
SparseMatrix interface to PETSc Mat.
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1521
virtual void swap(NumericVector< T > &v) libmesh_override
PetscErrorCode __libmesh_petsc_diff_solver_residual(SNES, Vec x, Vec r, void *ctx)
OStreamProxy out(std::cout)
std::unique_ptr< LinearSolutionMonitor > linear_solution_monitor
Definition: diff_solver.h:289
PetscErrorCode __libmesh_petsc_diff_solver_monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
Used for solving implicit systems of equations.