Actual source code: ex5.c


  2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
  3: This example tests PCVPBJacobiSetBlocks().\n\n";

  5: /*
  6:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
  7:    file automatically includes:
  8:      petscsys.h       - base PETSc routines   petscvec.h - vectors
  9:      petscmat.h - matrices
 10:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 11:      petscviewer.h - viewers               petscpc.h  - preconditioners
 12:      petscksp.h   - linear solvers
 13: */

 15: #include <petscsnes.h>

 17: /*
 18:    User-defined routines
 19: */
 20: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 21: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 22: extern PetscErrorCode FormInitialGuess(Vec);

 24: int main(int argc, char **argv)
 25: {
 26:   SNES        snes;       /* SNES context */
 27:   Vec         x, r, F, U; /* vectors */
 28:   Mat         J;          /* Jacobian matrix */
 29:   PetscInt    its, n = 5, i, maxit, maxf, lens[3] = {1, 2, 2};
 30:   PetscMPIInt size;
 31:   PetscScalar h, xp, v, none = -1.0;
 32:   PetscReal   abstol, rtol, stol, norm;
 33:   KSP         ksp;
 34:   PC          pc;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 38:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 39:   PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 40:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 41:   h = 1.0 / (n - 1);

 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44:      Create nonlinear solver context
 45:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 47:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 48:   PetscCall(SNESGetKSP(snes, &ksp));
 49:   PetscCall(KSPGetPC(ksp, &pc));
 50:   PetscCall(PCSetType(pc, PCVPBJACOBI));
 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Create vector data structures; set function evaluation routine
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   /*
 55:      Note that we form 1 vector from scratch and then duplicate as needed.
 56:   */
 57:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 58:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 59:   PetscCall(VecSetFromOptions(x));
 60:   PetscCall(VecDuplicate(x, &r));
 61:   PetscCall(VecDuplicate(x, &F));
 62:   PetscCall(VecDuplicate(x, &U));

 64:   /*
 65:      Set function evaluation routine and vector
 66:   */
 67:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));

 69:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 70:      Create matrix data structure; set Jacobian evaluation routine
 71:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 73:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
 74:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
 75:   PetscCall(MatSetFromOptions(J));
 76:   PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
 77:   PetscCall(MatSetVariableBlockSizes(J, 3, lens));

 79:   /*
 80:      Set Jacobian matrix data structure and default Jacobian evaluation
 81:      routine. User can override with:
 82:      -snes_fd : default finite differencing approximation of Jacobian
 83:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 84:                 (unless user explicitly sets preconditioner)
 85:      -snes_mf_operator : form preconditioning matrix as set by the user,
 86:                          but use matrix-free approx for Jacobian-vector
 87:                          products within Newton-Krylov method
 88:   */

 90:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Customize nonlinear solver; set runtime options
 94:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 96:   /*
 97:      Set names for some vectors to facilitate monitoring (optional)
 98:   */
 99:   PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
100:   PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));

102:   /*
103:      Set SNES/KSP/KSP/PC runtime options, e.g.,
104:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
105:   */
106:   PetscCall(SNESSetFromOptions(snes));

108:   /*
109:      Print parameters used for convergence testing (optional) ... just
110:      to demonstrate this routine; this information is also printed with
111:      the option -snes_view
112:   */
113:   PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
114:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Initialize application:
118:      Store right-hand-side of PDE and exact solution
119:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   xp = 0.0;
122:   for (i = 0; i < n; i++) {
123:     v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
124:     PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
125:     v = xp * xp * xp;
126:     PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
127:     xp += h;
128:   }

130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Evaluate initial guess; then solve nonlinear system
132:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133:   /*
134:      Note: The user should initialize the vector, x, with the initial guess
135:      for the nonlinear solver prior to calling SNESSolve().  In particular,
136:      to employ an initial guess of zero, the user should explicitly set
137:      this vector to zero by calling VecSet().
138:   */
139:   PetscCall(FormInitialGuess(x));
140:   PetscCall(SNESSolve(snes, NULL, x));
141:   PetscCall(SNESGetIterationNumber(snes, &its));
142:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:      Check solution and clean up
146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

148:   /*
149:      Check the error
150:   */
151:   PetscCall(VecAXPY(x, none, U));
152:   PetscCall(VecNorm(x, NORM_2, &norm));
153:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

155:   /*
156:      Free work space.  All PETSc objects should be destroyed when they
157:      are no longer needed.
158:   */
159:   PetscCall(VecDestroy(&x));
160:   PetscCall(VecDestroy(&r));
161:   PetscCall(VecDestroy(&U));
162:   PetscCall(VecDestroy(&F));
163:   PetscCall(MatDestroy(&J));
164:   PetscCall(SNESDestroy(&snes));
165:   PetscCall(PetscFinalize());
166:   return 0;
167: }

169: /*
170:    FormInitialGuess - Computes initial guess.

172:    Input/Output Parameter:
173: .  x - the solution vector
174: */
175: PetscErrorCode FormInitialGuess(Vec x)
176: {
177:   PetscScalar pfive = .50;

179:   PetscFunctionBeginUser;
180:   PetscCall(VecSet(x, pfive));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: /*
185:    FormFunction - Evaluates nonlinear function, F(x).

187:    Input Parameters:
188: .  snes - the SNES context
189: .  x - input vector
190: .  ctx - optional user-defined context, as set by SNESSetFunction()

192:    Output Parameter:
193: .  f - function vector

195:    Note:
196:    The user-defined context can contain any application-specific data
197:    needed for the function evaluation (such as various parameters, work
198:    vectors, and grid information).  In this program the context is just
199:    a vector containing the right-hand-side of the discretized PDE.
200:  */

202: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
203: {
204:   Vec                g = (Vec)ctx;
205:   const PetscScalar *xx, *gg;
206:   PetscScalar       *ff, d;
207:   PetscInt           i, n;

209:   PetscFunctionBeginUser;
210:   /*
211:      Get pointers to vector data.
212:        - For default PETSc vectors, VecGetArray() returns a pointer to
213:          the data array.  Otherwise, the routine is implementation dependent.
214:        - You MUST call VecRestoreArray() when you no longer need access to
215:          the array.
216:   */
217:   PetscCall(VecGetArrayRead(x, &xx));
218:   PetscCall(VecGetArray(f, &ff));
219:   PetscCall(VecGetArrayRead(g, &gg));

221:   /*
222:      Compute function
223:   */
224:   PetscCall(VecGetSize(x, &n));
225:   d     = (PetscReal)(n - 1);
226:   d     = d * d;
227:   ff[0] = xx[0];
228:   for (i = 1; i < n - 1; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - gg[i];
229:   ff[n - 1] = xx[n - 1] - 1.0;

231:   /*
232:      Restore vectors
233:   */
234:   PetscCall(VecRestoreArrayRead(x, &xx));
235:   PetscCall(VecRestoreArray(f, &ff));
236:   PetscCall(VecRestoreArrayRead(g, &gg));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: /*
241:    FormJacobian - Evaluates Jacobian matrix.

243:    Input Parameters:
244: .  snes - the SNES context
245: .  x - input vector
246: .  dummy - optional user-defined context (not used here)

248:    Output Parameters:
249: .  jac - Jacobian matrix
250: .  B - optionally different preconditioning matrix

252: */

254: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
255: {
256:   const PetscScalar *xx;
257:   PetscScalar        A[3], d;
258:   PetscInt           i, n, j[3];

260:   PetscFunctionBeginUser;
261:   /*
262:      Get pointer to vector data
263:   */
264:   PetscCall(VecGetArrayRead(x, &xx));

266:   /*
267:      Compute Jacobian entries and insert into matrix.
268:       - Note that in this case we set all elements for a particular
269:         row at once.
270:   */
271:   PetscCall(VecGetSize(x, &n));
272:   d = (PetscReal)(n - 1);
273:   d = d * d;

275:   /*
276:      Interior grid points
277:   */
278:   for (i = 1; i < n - 1; i++) {
279:     j[0] = i - 1;
280:     j[1] = i;
281:     j[2] = i + 1;
282:     A[0] = d;
283:     A[1] = -2.0 * d + 2.0 * xx[i];
284:     A[2] = d;
285:     PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
286:   }

288:   /*
289:      Boundary points
290:   */
291:   i    = 0;
292:   A[0] = 1.0;

294:   PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));

296:   i    = n - 1;
297:   A[0] = 1.0;

299:   PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));

301:   /*
302:      Restore vector
303:   */
304:   PetscCall(VecRestoreArrayRead(x, &xx));

306:   /*
307:      Assemble matrix
308:   */
309:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
310:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
311:   if (jac != B) {
312:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
313:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
314:   }
315:   PetscFunctionReturn(PETSC_SUCCESS);
316: }

318: /*TEST

320:    testset:
321:      args: -snes_monitor_short -snes_view -ksp_monitor
322:      output_file: output/ex5_1.out
323:      filter: grep -v "type: seqaij"

325:      test:
326:       suffix: 1

328:      test:
329:       suffix: cuda
330:       requires: cuda
331:       args: -mat_type aijcusparse -vec_type cuda

333:      test:
334:       suffix: kok
335:       requires: kokkos_kernels
336:       args: -mat_type aijkokkos -vec_type kokkos

338:    # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly
339:    # the solution is wrong on purpose
340:    test:
341:       requires: !single !complex
342:       suffix: transpose_only
343:       args: -snes_monitor_short -snes_view -ksp_monitor -snes_type ksptransposeonly -pc_type ilu -snes_test_jacobian -snes_test_jacobian_view -ksp_view_rhs -ksp_view_solution -ksp_view_mat_explicit -ksp_view_preconditioned_operator_explicit

345: TEST*/