Actual source code: ex78.c


  2: static char help[] = "Newton methods to solve u''  = f in parallel with periodic boundary conditions.\n\n";

  4: /*
  5:    Compare this example to ex3.c that handles Dirichlet boundary conditions

  7:    Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate
  8:    handling periodic boundary conditions for nonlinear problems.

 10:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 11:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 12:    file automatically includes:
 13:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 14:      petscmat.h - matrices
 15:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 16:      petscviewer.h - viewers               petscpc.h  - preconditioners
 17:      petscksp.h   - linear solvers
 18: */

 20: #include <petscdm.h>
 21: #include <petscdmda.h>
 22: #include <petscsnes.h>

 24: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 25: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);

 27: int main(int argc, char **argv)
 28: {
 29:   SNES         snes; /* SNES context */
 30:   Mat          J;    /* Jacobian matrix */
 31:   DM           da;
 32:   Vec          x, r; /* vectors */
 33:   PetscInt     N = 5;
 34:   MatNullSpace constants;

 36:   PetscFunctionBeginUser;
 37:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 38:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));

 40:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41:      Create nonlinear solver context
 42:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 44:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 47:      Create vector data structures; set function evaluation routine
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   /*
 51:      Create distributed array (DMDA) to manage parallel grid and vectors
 52:   */
 53:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, N, 1, 1, NULL, &da));
 54:   PetscCall(DMSetFromOptions(da));
 55:   PetscCall(DMSetUp(da));

 57:   /*
 58:      Extract global and local vectors from DMDA; then duplicate for remaining
 59:      vectors that are the same types
 60:   */
 61:   PetscCall(DMCreateGlobalVector(da, &x));
 62:   PetscCall(VecDuplicate(x, &r));

 64:   /*
 65:      Set function evaluation routine and vector.  Whenever the nonlinear
 66:      solver needs to compute the nonlinear function, it will call this
 67:      routine.
 68:       - Note that the final routine argument is the user-defined
 69:         context that provides application-specific data for the
 70:         function evaluation routine.
 71:   */
 72:   PetscCall(SNESSetFunction(snes, r, FormFunction, da));

 74:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75:      Create matrix data structure; set Jacobian evaluation routine
 76:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 77:   PetscCall(DMCreateMatrix(da, &J));
 78:   PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &constants));
 79:   PetscCall(MatSetNullSpace(J, constants));
 80:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, da));

 82:   PetscCall(SNESSetFromOptions(snes));
 83:   PetscCall(SNESSolve(snes, NULL, x));

 85:   PetscCall(VecDestroy(&x));
 86:   PetscCall(VecDestroy(&r));
 87:   PetscCall(MatDestroy(&J));
 88:   PetscCall(MatNullSpaceDestroy(&constants));
 89:   PetscCall(SNESDestroy(&snes));
 90:   PetscCall(DMDestroy(&da));
 91:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }

 95: /*
 96:    FormFunction - Evaluates nonlinear function, F(x).

 98:    Input Parameters:
 99: .  snes - the SNES context
100: .  x - input vector
101: .  ctx - optional user-defined context, as set by SNESSetFunction()

103:    Output Parameter:
104: .  f - function vector

106:    Note:
107:    The user-defined context can contain any application-specific
108:    data needed for the function evaluation.
109: */
110: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
111: {
112:   DM           da = (DM)ctx;
113:   PetscScalar *xx, *ff;
114:   PetscReal    h;
115:   PetscInt     i, M, xs, xm;
116:   Vec          xlocal;

118:   PetscFunctionBeginUser;
119:   /* Get local work vector */
120:   PetscCall(DMGetLocalVector(da, &xlocal));

122:   /*
123:      Scatter ghost points to local vector, using the 2-step process
124:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
125:      By placing code between these two statements, computations can
126:      be done while messages are in transition.
127:   */
128:   PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal));
129:   PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal));

131:   /*
132:      Get pointers to vector data.
133:        - The vector xlocal includes ghost point; the vectors x and f do
134:          NOT include ghost points.
135:        - Using DMDAVecGetArray() allows accessing the values using global ordering
136:   */
137:   PetscCall(DMDAVecGetArray(da, xlocal, &xx));
138:   PetscCall(DMDAVecGetArray(da, f, &ff));

140:   /*
141:      Get local grid boundaries (for 1-dimensional DMDA):
142:        xs, xm  - starting grid index, width of local grid (no ghost points)
143:   */
144:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
145:   PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));

147:   /*
148:      Compute function over locally owned part of the grid
149:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
150:   */
151:   h = 1.0 / M;
152:   for (i = xs; i < xs + xm; i++) ff[i] = (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) / (h * h) - PetscSinReal(2.0 * PETSC_PI * i * h);

154:   /*
155:      Restore vectors
156:   */
157:   PetscCall(DMDAVecRestoreArray(da, xlocal, &xx));
158:   PetscCall(DMDAVecRestoreArray(da, f, &ff));
159:   PetscCall(DMRestoreLocalVector(da, &xlocal));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }
162: /* ------------------------------------------------------------------- */
163: /*
164:    FormJacobian - Evaluates Jacobian matrix.

166:    Input Parameters:
167: .  snes - the SNES context
168: .  x - input vector
169: .  dummy - optional user-defined context (not used here)

171:    Output Parameters:
172: .  jac - Jacobian matrix
173: .  B - optionally different preconditioning matrix
174: .  flag - flag indicating matrix structure
175: */
176: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *ctx)
177: {
178:   PetscScalar *xx, A[3];
179:   PetscInt     i, M, xs, xm;
180:   DM           da = (DM)ctx;
181:   MatStencil   row, cols[3];
182:   PetscReal    h;

184:   PetscFunctionBeginUser;
185:   /*
186:      Get pointer to vector data
187:   */
188:   PetscCall(DMDAVecGetArrayRead(da, x, &xx));
189:   PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

191:   /*
192:     Get range of locally owned matrix
193:   */
194:   PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));

196:   PetscCall(MatZeroEntries(jac));
197:   h = 1.0 / M;
198:   /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
199:   for (i = xs; i < xs + xm; i++) {
200:     row.i     = i;
201:     cols[0].i = i - 1;
202:     cols[1].i = i;
203:     cols[2].i = i + 1;
204:     A[0] = A[2] = 1.0 / (h * h);
205:     A[1]        = -2.0 / (h * h);
206:     PetscCall(MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES));
207:   }

209:   PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
210:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
211:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }

215: /*TEST

217:    test:
218:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
219:       requires: !single

221:    test:
222:       suffix: 2
223:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc
224:       requires: !single

226:    test:
227:       suffix: 3
228:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false
229:       requires: !single

231: TEST*/