petsc_diff_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
25 
26 #ifdef LIBMESH_HAVE_PETSC
27 
28 namespace libMesh
29 {
30 
31 //--------------------------------------------------------------------
32 // Functions with C linkage to pass to PETSc. PETSc will call these
33 // methods as needed.
34 //
35 // Since they must have C linkage they have no knowledge of a namespace.
36 // Give them an obscure name to avoid namespace pollution.
37 extern "C"
38 {
39  // Older versions of PETSc do not have the different int typedefs.
40  // On 64-bit machines, PetscInt may actually be a long long int.
41  // This change occurred in Petsc-2.2.1.
42 #if PETSC_VERSION_LESS_THAN(2,2,1)
43  typedef int PetscErrorCode;
44  typedef int PetscInt;
45 #endif
46 
47 // Function to hand to PETSc's SNES,
48 // which monitors convergence at X
49 PetscErrorCode
50 __libmesh_petsc_diff_solver_monitor (SNES snes, PetscInt its,
51  PetscReal fnorm, void *ctx)
52 {
53  PetscDiffSolver& solver =
54  *(static_cast<PetscDiffSolver*> (ctx));
55 
56  if (solver.verbose)
57  libMesh::out << " PetscDiffSolver step " << its
58  << ", |residual|_2 = " << fnorm << std::endl;
59  if (solver.linear_solution_monitor.get()) {
60  int ierr = 0;
61 
62  Vec petsc_delta_u;
63  ierr = SNESGetSolutionUpdate(snes, &petsc_delta_u);
65  PetscVector<Number> delta_u(petsc_delta_u, solver.comm());
66  delta_u.close();
67 
68  Vec petsc_u;
69  ierr = SNESGetSolution(snes, &petsc_u);
71  PetscVector<Number> u(petsc_u, solver.comm());
72  u.close();
73 
74  Vec petsc_res;
75  ierr = SNESGetFunction(snes, &petsc_res, NULL, NULL);
77  PetscVector<Number> res(petsc_res, solver.comm());
78  res.close();
79 
80  (*solver.linear_solution_monitor)(
81  delta_u, delta_u.l2_norm(),
82  u, u.l2_norm(),
83  res, res.l2_norm(), its);
84  }
85  return 0;
86 }
87 
88 // Functions to hand to PETSc's SNES,
89 // which compute the residual or jacobian at X
90 PetscErrorCode
91 __libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void *ctx)
92 {
93  libmesh_assert(x);
94  libmesh_assert(r);
95  libmesh_assert(ctx);
96 
97  PetscDiffSolver& solver =
98  *(static_cast<PetscDiffSolver*> (ctx));
99  ImplicitSystem &sys = solver.system();
100 
101  if (solver.verbose)
102  libMesh::out << "Assembling the residual" << std::endl;
103 
104  PetscVector<Number>& X_system =
105  *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
106  PetscVector<Number>& R_system =
107  *libmesh_cast_ptr<PetscVector<Number>*>(sys.rhs);
108  PetscVector<Number> X_input(x, sys.comm()), R_input(r, sys.comm());
109 
110  // DiffSystem assembles from the solution and into the rhs, so swap
111  // those with our input vectors before assembling. They'll probably
112  // already be references to the same vectors, but PETSc might do
113  // something tricky.
114  X_input.swap(X_system);
115  R_input.swap(R_system);
116 
117  // We may need to correct a non-conforming solution
119 
120  // We may need to localize a parallel solution
121  sys.update();
122 
123  // Do DiffSystem assembly
124  sys.assembly(true, false);
125  R_system.close();
126 
127  // Swap back
128  X_input.swap(X_system);
129  R_input.swap(R_system);
130 
131  // No errors, we hope
132  return 0;
133 }
134 
135 
136 PetscErrorCode
137 __libmesh_petsc_diff_solver_jacobian (SNES, Vec x, Mat *libmesh_dbg_var(j), Mat *pc,
138  MatStructure *msflag, void *ctx)
139 {
140  libmesh_assert(x);
141  libmesh_assert(j);
142 // libmesh_assert_equal_to (pc, j); // We don't use separate preconditioners yet
143  libmesh_assert(ctx);
144 
145  PetscDiffSolver& solver =
146  *(static_cast<PetscDiffSolver*> (ctx));
147  ImplicitSystem &sys = solver.system();
148 
149  if (solver.verbose)
150  libMesh::out << "Assembling the Jacobian" << std::endl;
151 
152  PetscVector<Number>& X_system =
153  *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
154  PetscVector<Number> X_input(x, sys.comm());
155 
156  PetscMatrix<Number> J_input(*pc, sys.comm());
157  PetscMatrix<Number>& J_system =
158  *libmesh_cast_ptr<PetscMatrix<Number>*>(sys.matrix);
159 
160  // DiffSystem assembles from the solution and into the jacobian, so
161  // swap those with our input vectors before assembling. They'll
162  // probably already be references to the same vectors, but PETSc
163  // might do something tricky.
164  X_input.swap(X_system);
165  J_input.swap(J_system);
166 
167  // We may need to correct a non-conforming solution
169 
170  // We may need to localize a parallel solution
171  sys.update();
172 
173  // Do DiffSystem assembly
174  sys.assembly(false, true);
175  J_system.close();
176 
177  // Swap back
178  X_input.swap(X_system);
179  J_input.swap(J_system);
180 
181  *msflag = SAME_NONZERO_PATTERN;
182 
183  // No errors, we hope
184  return 0;
185 }
186 
187 } // extern "C"
188 
189 
191  : Parent(s)
192 {
193 }
194 
195 
197 {
198  START_LOG("init()", "PetscDiffSolver");
199 
200  Parent::init();
201 
202  int ierr=0;
203 
204 #if PETSC_VERSION_LESS_THAN(2,1,2)
205  // At least until Petsc 2.1.1, the SNESCreate had a different
206  // calling syntax. The second argument was of type SNESProblemType,
207  // and could have a value of either SNES_NONLINEAR_EQUATIONS or
208  // SNES_UNCONSTRAINED_MINIMIZATION.
209  ierr = SNESCreate(this->comm().get(), SNES_NONLINEAR_EQUATIONS, &_snes);
210  LIBMESH_CHKERRABORT(ierr);
211 #else
212  ierr = SNESCreate(this->comm().get(),&_snes);
213  LIBMESH_CHKERRABORT(ierr);
214 #endif
215 
216 #if PETSC_VERSION_LESS_THAN(2,3,3)
217  ierr = SNESSetMonitor (_snes, __libmesh_petsc_diff_solver_monitor,
218  this, PETSC_NULL);
219 #else
220  // API name change in PETSc 2.3.3
221  ierr = SNESMonitorSet (_snes, __libmesh_petsc_diff_solver_monitor,
222  this, PETSC_NULL);
223 #endif
224  LIBMESH_CHKERRABORT(ierr);
225 
226  ierr = SNESSetFromOptions(_snes);
227  LIBMESH_CHKERRABORT(ierr);
228 
229  STOP_LOG("init()", "PetscDiffSolver");
230 }
231 
232 
233 
235 {
236 }
237 
238 
239 
241 {
242  START_LOG("clear()", "PetscDiffSolver");
243 
244  int ierr=0;
245 
246  ierr = LibMeshSNESDestroy(&_snes);
247  LIBMESH_CHKERRABORT(ierr);
248 
249  STOP_LOG("clear()", "PetscDiffSolver");
250 }
251 
252 
253 
255 {
256  Parent::reinit();
257 }
258 
259 
260 
262 {
263  switch (r)
264  {
265  case SNES_CONVERGED_FNORM_ABS:
267  case SNES_CONVERGED_FNORM_RELATIVE:
269 #if PETSC_VERSION_LESS_THAN(3,2,1)
270  case SNES_CONVERGED_PNORM_RELATIVE:
271 #else
272  case SNES_CONVERGED_SNORM_RELATIVE:
273 #endif
275 #if !PETSC_VERSION_LESS_THAN(2,3,3)
276  case SNES_CONVERGED_ITS:
277 #endif
278  case SNES_CONVERGED_TR_DELTA:
280  case SNES_DIVERGED_FUNCTION_DOMAIN:
281  case SNES_DIVERGED_FUNCTION_COUNT:
282  case SNES_DIVERGED_FNORM_NAN:
283 #if !PETSC_VERSION_LESS_THAN(3,3,0)
284  case SNES_DIVERGED_INNER:
285 #endif
286 #if !PETSC_VERSION_LESS_THAN(2,3,2)
287  case SNES_DIVERGED_LINEAR_SOLVE:
288 #endif
289  case SNES_DIVERGED_LOCAL_MIN:
291  case SNES_DIVERGED_MAX_IT:
293 #if !PETSC_VERSION_LESS_THAN(3,2,0)
294  case SNES_DIVERGED_LINE_SEARCH:
296 #endif
297  // In PETSc, SNES_CONVERGED_ITERATING means
298  // the solve is still iterating, but by the
299  // time we get here, we must have either
300  // converged or diverged, so
301  // SNES_CONVERGED_ITERATING is invalid.
302  case SNES_CONVERGED_ITERATING:
304  }
306 }
307 
308 
309 
311 {
312  this->init();
313 
314  START_LOG("solve()", "PetscDiffSolver");
315 
317  *(libmesh_cast_ptr<PetscVector<Number>*>(_system.solution.get()));
318  PetscMatrix<Number> &jac =
319  *(libmesh_cast_ptr<PetscMatrix<Number>*>(_system.matrix));
321  *(libmesh_cast_ptr<PetscVector<Number>*>(_system.rhs));
322 
323 #ifdef LIBMESH_ENABLE_CONSTRAINTS
325 #endif
326 
327  int ierr = 0;
328 
329  ierr = SNESSetFunction (_snes, r.vec(),
331  LIBMESH_CHKERRABORT(ierr);
332 
333  ierr = SNESSetJacobian (_snes, jac.mat(), jac.mat(),
335  LIBMESH_CHKERRABORT(ierr);
336 
337 # if PETSC_VERSION_LESS_THAN(2,2,0)
338 
339  ierr = SNESSolve (_snes, x.vec(), &_outer_iterations);
340  LIBMESH_CHKERRABORT(ierr);
341 
342 // 2.2.x style
343 #elif PETSC_VERSION_LESS_THAN(2,3,0)
344 
345  ierr = SNESSolve (_snes, x.vec());
346  LIBMESH_CHKERRABORT(ierr);
347 
348 // 2.3.x & newer style
349 #else
350 
351  ierr = SNESSolve (_snes, PETSC_NULL, x.vec());
352  LIBMESH_CHKERRABORT(ierr);
353 
354 #endif
355 
356  STOP_LOG("solve()", "PetscDiffSolver");
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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo