Actual source code: ex5.c

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

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

 14: #include <petscsnes.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

107:   /*
108:      Print parameters used for convergence testing (optional) ... just
109:      to demonstrate this routine; this information is also printed with
110:      the option -snes_view
111:   */
112:   PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
113:   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));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

220:   /*
221:      Compute function
222:   */
223:   PetscCall(VecGetSize(x, &n));
224:   d     = (PetscReal)(n - 1);
225:   d     = d * d;
226:   ff[0] = xx[0];
227:   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];
228:   ff[n - 1] = xx[n - 1] - 1.0;

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

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

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

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

251: */

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

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

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

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

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

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

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

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

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

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

317: /*TEST

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

324:      test:
325:       suffix: 1

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

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

337:    # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly
338:    # the solution is wrong on purpose
339:    test:
340:       requires: !single !complex
341:       suffix: transpose_only
342:       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

344: TEST*/