Actual source code: ex78.c


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

  4: /*T
  5:    Concepts: SNES^basic parallel example
  6:    Concepts: periodic boundary conditions
  7:    Processors: n
  8: T*/

 10: /*
 11:    Compare this example to ex3.c that handles Dirichlet boundary conditions

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

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

 26: #include <petscdm.h>
 27: #include <petscdmda.h>
 28: #include <petscsnes.h>

 30: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 31: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);

 33: int main(int argc,char **argv)
 34: {
 35:   SNES           snes;                 /* SNES context */
 36:   Mat            J;                    /* Jacobian matrix */
 37:   DM             da;
 38:   Vec            x,r;              /* vectors */
 39:   PetscInt       N = 5;
 40:   MatNullSpace   constants;

 42:   PetscInitialize(&argc,&argv,(char*)0,help);
 43:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);

 45:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 46:      Create nonlinear solver context
 47:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 49:   SNESCreate(PETSC_COMM_WORLD,&snes);

 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Create vector data structures; set function evaluation routine
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 55:   /*
 56:      Create distributed array (DMDA) to manage parallel grid and vectors
 57:   */
 58:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);
 59:   DMSetFromOptions(da);
 60:   DMSetUp(da);

 62:   /*
 63:      Extract global and local vectors from DMDA; then duplicate for remaining
 64:      vectors that are the same types
 65:   */
 66:   DMCreateGlobalVector(da,&x);
 67:   VecDuplicate(x,&r);

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

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Create matrix data structure; set Jacobian evaluation routine
 81:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 82:   DMCreateMatrix(da,&J);
 83:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
 84:   MatSetNullSpace(J,constants);
 85:   SNESSetJacobian(snes,J,J,FormJacobian,da);

 87:   SNESSetFromOptions(snes);
 88:   SNESSolve(snes,NULL,x);

 90:   VecDestroy(&x);
 91:   VecDestroy(&r);
 92:   MatDestroy(&J);
 93:   MatNullSpaceDestroy(&constants);
 94:   SNESDestroy(&snes);
 95:   DMDestroy(&da);
 96:   PetscFinalize();
 97:   return 0;
 98: }

100: /*
101:    FormFunction - Evaluates nonlinear function, F(x).

103:    Input Parameters:
104: .  snes - the SNES context
105: .  x - input vector
106: .  ctx - optional user-defined context, as set by SNESSetFunction()

108:    Output Parameter:
109: .  f - function vector

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

124:   /* Get local work vector */
125:   DMGetLocalVector(da,&xlocal);

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

136:   /*
137:      Get pointers to vector data.
138:        - The vector xlocal includes ghost point; the vectors x and f do
139:          NOT include ghost points.
140:        - Using DMDAVecGetArray() allows accessing the values using global ordering
141:   */
142:   DMDAVecGetArray(da,xlocal,&xx);
143:   DMDAVecGetArray(da,f,&ff);

145:   /*
146:      Get local grid boundaries (for 1-dimensional DMDA):
147:        xs, xm  - starting grid index, width of local grid (no ghost points)
148:   */
149:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
150:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

152:   /*
153:      Compute function over locally owned part of the grid
154:      Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
155:   */
156:   h = 1.0/M;
157:   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);

159:   /*
160:      Restore vectors
161:   */
162:   DMDAVecRestoreArray(da,xlocal,&xx);
163:   DMDAVecRestoreArray(da,f,&ff);
164:   DMRestoreLocalVector(da,&xlocal);
165:   return 0;
166: }
167: /* ------------------------------------------------------------------- */
168: /*
169:    FormJacobian - Evaluates Jacobian matrix.

171:    Input Parameters:
172: .  snes - the SNES context
173: .  x - input vector
174: .  dummy - optional user-defined context (not used here)

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

190:   /*
191:      Get pointer to vector data
192:   */
193:   DMDAVecGetArrayRead(da,x,&xx);
194:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

196:   /*
197:     Get range of locally owned matrix
198:   */
199:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

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

213:   DMDAVecRestoreArrayRead(da,x,&xx);
214:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
215:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
216:   return 0;
217: }

219: /*TEST

221:    test:
222:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
223:       requires: !single

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

230:    test:
231:       suffix: 3
232:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false
233:       requires: !single

235: TEST*/