Actual source code: ex78.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

146:   /*
147:      Compute function over locally owned part of the grid
148:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
149:   */
150:   h = 1.0 / M;
151:   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);

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

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

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

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

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

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

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

214: /*TEST

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

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

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

230: TEST*/