petsc_diff_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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  {
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 
93  *(static_cast<PetscDiffSolver*> (ctx));
94  ImplicitSystem & sys = solver.system();
95 
96  if (solver.verbose)
97  libMesh::out << "Assembling the residual" << std::endl;
98 
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 correct a non-conforming solution
114 
115  // We may need to localize a parallel solution
116  sys.update();
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 #if PETSC_RELEASE_LESS_THAN(3,5,0)
132  PetscErrorCode
134  Vec x,
135  Mat * libmesh_dbg_var(j),
136  Mat * pc,
137  MatStructure * msflag,
138  void * ctx)
139 #else
140  PetscErrorCode
142  Vec x,
143  Mat libmesh_dbg_var(j),
144  Mat pc,
145  void * ctx)
146 #endif
147  {
148  libmesh_assert(x);
149  libmesh_assert(j);
150  // libmesh_assert_equal_to (pc, j); // We don't use separate preconditioners yet
151  libmesh_assert(ctx);
152 
154  *(static_cast<PetscDiffSolver*> (ctx));
155  ImplicitSystem & sys = solver.system();
156 
157  if (solver.verbose)
158  libMesh::out << "Assembling the Jacobian" << std::endl;
159 
161  *cast_ptr<PetscVector<Number> *>(sys.solution.get());
162  PetscVector<Number> X_input(x, sys.comm());
163 
164 #if PETSC_RELEASE_LESS_THAN(3,5,0)
165  PetscMatrix<Number> J_input(*pc, sys.comm());
166 #else
167  PetscMatrix<Number> J_input(pc, sys.comm());
168 #endif
170  *cast_ptr<PetscMatrix<Number> *>(sys.matrix);
171 
172  // DiffSystem assembles from the solution and into the jacobian, so
173  // swap those with our input vectors before assembling. They'll
174  // probably already be references to the same vectors, but PETSc
175  // might do something tricky.
176  X_input.swap(X_system);
177  J_input.swap(J_system);
178 
179  // We may need to correct a non-conforming solution
181 
182  // We may need to localize a parallel solution
183  sys.update();
184 
185  // Do DiffSystem assembly
186  sys.assembly(false, true);
187  J_system.close();
188 
189  // Swap back
190  X_input.swap(X_system);
191  J_input.swap(J_system);
192 
193 #if PETSC_RELEASE_LESS_THAN(3,5,0)
194  *msflag = SAME_NONZERO_PATTERN;
195 #endif
196  // No errors, we hope
197  return 0;
198  }
199 
200 } // extern "C"
201 
202 
204  : Parent(s)
205 {
206 }
207 
208 
210 {
211  LOG_SCOPE("init()", "PetscDiffSolver");
212 
213  Parent::init();
214 
215  int ierr=0;
216 
217  ierr = SNESCreate(this->comm().get(),&_snes);
218  LIBMESH_CHKERR(ierr);
219 
220  ierr = SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor,
221  this, PETSC_NULL);
222  LIBMESH_CHKERR(ierr);
223 
224  if (libMesh::on_command_line("--solver_system_names"))
225  {
226  ierr = SNESSetOptionsPrefix(_snes, (_system.name()+"_").c_str());
227  LIBMESH_CHKERR(ierr);
228  }
229 
230  ierr = SNESSetFromOptions(_snes);
231  LIBMESH_CHKERR(ierr);
232 
233  KSP my_ksp;
234  ierr = SNESGetKSP(_snes, &my_ksp);
235  LIBMESH_CHKERR(ierr);
236 
237  PC my_pc;
238  ierr = KSPGetPC(my_ksp, &my_pc);
239  LIBMESH_CHKERR(ierr);
240 
242 }
243 
244 
245 
247 {
248 }
249 
250 
251 
253 {
254  LOG_SCOPE("clear()", "PetscDiffSolver");
255 
256  int ierr = LibMeshSNESDestroy(&_snes);
257  LIBMESH_CHKERR(ierr);
258 }
259 
260 
261 
263 {
264  Parent::reinit();
265 
266  KSP my_ksp;
267  int ierr = SNESGetKSP(_snes, &my_ksp);
268  LIBMESH_CHKERR(ierr);
269 
270  PC my_pc;
271  ierr = KSPGetPC(my_ksp, &my_pc);
272  LIBMESH_CHKERR(ierr);
273 
275 }
276 
277 
278 
280 {
281  switch (r)
282  {
283  case SNES_CONVERGED_FNORM_ABS:
285  case SNES_CONVERGED_FNORM_RELATIVE:
287 #if PETSC_VERSION_LESS_THAN(3,2,1)
288  case SNES_CONVERGED_PNORM_RELATIVE:
289 #else
290  case SNES_CONVERGED_SNORM_RELATIVE:
291 #endif
293  case SNES_CONVERGED_ITS:
294  case SNES_CONVERGED_TR_DELTA:
296  case SNES_DIVERGED_FUNCTION_DOMAIN:
297  case SNES_DIVERGED_FUNCTION_COUNT:
298  case SNES_DIVERGED_FNORM_NAN:
299 #if !PETSC_VERSION_LESS_THAN(3,3,0)
300  case SNES_DIVERGED_INNER:
301 #endif
302  case SNES_DIVERGED_LINEAR_SOLVE:
303  case SNES_DIVERGED_LOCAL_MIN:
305  case SNES_DIVERGED_MAX_IT:
307 #if PETSC_VERSION_LESS_THAN(3,2,0)
308  case SNES_DIVERGED_LS_FAILURE:
309 #else
310  case SNES_DIVERGED_LINE_SEARCH:
311 #endif
313  // In PETSc, SNES_CONVERGED_ITERATING means
314  // the solve is still iterating, but by the
315  // time we get here, we must have either
316  // converged or diverged, so
317  // SNES_CONVERGED_ITERATING is invalid.
318  case SNES_CONVERGED_ITERATING:
320  default:
321  break;
322  }
324 }
325 
326 
327 
329 {
330  this->init();
331 
332  LOG_SCOPE("solve()", "PetscDiffSolver");
333 
334  PetscVector<Number> & x =
335  *(cast_ptr<PetscVector<Number> *>(_system.solution.get()));
336  PetscMatrix<Number> & jac =
337  *(cast_ptr<PetscMatrix<Number> *>(_system.matrix));
338  PetscVector<Number> & r =
339  *(cast_ptr<PetscVector<Number> *>(_system.rhs));
340 
341 #ifdef LIBMESH_ENABLE_CONSTRAINTS
343 #endif
344 
345  int ierr = 0;
346 
347  ierr = SNESSetFunction (_snes, r.vec(),
349  LIBMESH_CHKERR(ierr);
350 
351  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
353  LIBMESH_CHKERR(ierr);
354 
355  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
356  LIBMESH_CHKERR(ierr);
357 
358  SNESConvergedReason reason;
359  SNESGetConvergedReason(_snes, &reason);
360 
361  this->clear();
362 
363  return convert_solve_result(reason);
364 }
365 
366 
367 } // namespace libMesh
368 
369 #endif // LIBMESH_HAVE_PETSC
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
PetscErrorCode __libmesh_petsc_diff_solver_jacobian(SNES, Vec x, Mat *libmesh_dbg_var(j), Mat *pc, MatStructure *msflag, void *ctx) PetscErrorCode __libmesh_petsc_diff_solver_jacobian(SNES
virtual void reinit() libmesh_override
ImplicitSystem & sys
virtual void init() libmesh_override
UniquePtr< NumericVector< Number > > current_local_solution
Definition: system.h:1533
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
PetscErrorCode Vec Mat Mat pc
const std::string & name() const
Definition: system.h:1996
PetscDiffSolver(sys_type &system)
libmesh_assert(j)
PetscDiffSolver & solver
PetscMatrix< Number > & J_system
virtual void init()
Definition: diff_solver.C:69
virtual void assembly(bool, bool, bool=false, bool=false)
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2028
sys_type & _system
Definition: diff_solver.h:318
virtual void close() libmesh_override
Definition: petsc_matrix.C:897
const sys_type & system() const
Definition: diff_solver.h:137
PetscErrorCode Vec Mat Mat void * ctx
UniquePtr< NumericVector< Number > > solution
Definition: system.h:1521
if(solver.verbose) libMesh PetscVector< Number > & X_system
PetscErrorCode Vec Mat libmesh_dbg_var(j)
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
PetscMatrix< Number > J_input(pc, sys.comm())
virtual void close() libmesh_override
Definition: petsc_vector.h:807
virtual void update()
Definition: system.C:423
virtual unsigned int solve() libmesh_override
PetscErrorCode ierr
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
const Parallel::Communicator & comm() const
SparseMatrix< Number > * matrix
SparseMatrix interface to PETSc Mat.
UniquePtr< LinearSolutionMonitor > linear_solution_monitor
Definition: diff_solver.h:288
PetscErrorCode __libmesh_petsc_diff_solver_residual(SNES, Vec x, Vec r, void *ctx)
OStreamProxy out(std::cout)
PetscVector< Number > X_input(x, sys.comm())
PetscErrorCode __libmesh_petsc_diff_solver_monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
Used for solving implicit systems of equations.