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;

 37:   PetscInitialize(&argc, &argv, (char *)0, help);
 38:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 40:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 41:   h = 1.0 / (n - 1);

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

 47:   SNESCreate(PETSC_COMM_WORLD, &snes);
 48:   SNESGetKSP(snes, &ksp);
 49:   KSPGetPC(ksp, &pc);
 50:   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:   VecCreate(PETSC_COMM_WORLD, &x);
 58:   VecSetSizes(x, PETSC_DECIDE, n);
 59:   VecSetFromOptions(x);
 60:   VecDuplicate(x, &r);
 61:   VecDuplicate(x, &F);
 62:   VecDuplicate(x, &U);

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

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

 73:   MatCreate(PETSC_COMM_WORLD, &J);
 74:   MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n);
 75:   MatSetFromOptions(J);
 76:   MatSeqAIJSetPreallocation(J, 3, NULL);
 77:   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:   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:   PetscObjectSetName((PetscObject)x, "Approximate Solution");
100:   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:   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:   SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf);
114:   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:     VecSetValues(F, 1, &i, &v, INSERT_VALUES);
125:     v = xp * xp * xp;
126:     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:   FormInitialGuess(x);
140:   SNESSolve(snes, NULL, x);
141:   SNESGetIterationNumber(snes, &its);
142:   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:   VecAXPY(x, none, U);
152:   VecNorm(x, NORM_2, &norm);
153:   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:   VecDestroy(&x);
160:   VecDestroy(&r);
161:   VecDestroy(&U);
162:   VecDestroy(&F);
163:   MatDestroy(&J);
164:   SNESDestroy(&snes);
165:   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;
178:   VecSet(x, pfive);
179:   return 0;
180: }

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

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

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

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

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

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

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

228:   /*
229:      Restore vectors
230:   */
231:   VecRestoreArrayRead(x, &xx);
232:   VecRestoreArray(f, &ff);
233:   VecRestoreArrayRead(g, &gg);
234:   return 0;
235: }

237: /*
238:    FormJacobian - Evaluates Jacobian matrix.

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

245:    Output Parameters:
246: .  jac - Jacobian matrix
247: .  B - optionally different preconditioning matrix

249: */

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

257:   /*
258:      Get pointer to vector data
259:   */
260:   VecGetArrayRead(x, &xx);

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

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

283:   /*
284:      Boundary points
285:   */
286:   i    = 0;
287:   A[0] = 1.0;

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

291:   i    = n - 1;
292:   A[0] = 1.0;

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

296:   /*
297:      Restore vector
298:   */
299:   VecRestoreArrayRead(x, &xx);

301:   /*
302:      Assemble matrix
303:   */
304:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
305:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
306:   if (jac != B) {
307:     MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
308:     MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
309:   }
310:   return 0;
311: }

313: /*TEST

315:    testset:
316:      args: -snes_monitor_short -snes_view -ksp_monitor
317:      output_file: output/ex5_1.out
318:      filter: grep -v "type: seqaij"

320:      test:
321:       suffix: 1

323:      test:
324:       suffix: cuda
325:       requires: cuda
326:       args: -mat_type aijcusparse -vec_type cuda

328:      test:
329:       suffix: kok
330:       requires: kokkos_kernels
331:       args: -mat_type aijkokkos -vec_type kokkos

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

340: TEST*/