Actual source code: ex78.c

petsc-3.9.3 2018-07-02
Report Typos and Errors

  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*/



 12: /*
 13:    Compare this example to ex3.c that handles Dirichlet boundary conditions

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

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

 28:  #include <petscdm.h>
 29:  #include <petscdmda.h>
 30:  #include <petscsnes.h>


 33: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 34: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);

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

 46:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 47:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create nonlinear solver context
 51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 53:   SNESCreate(PETSC_COMM_WORLD,&snes);

 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Create vector data structures; set function evaluation routine
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 59:   /*
 60:      Create distributed array (DMDA) to manage parallel grid and vectors
 61:   */
 62:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);
 63:   DMSetFromOptions(da);
 64:   DMSetUp(da);

 66:   /*
 67:      Extract global and local vectors from DMDA; then duplicate for remaining
 68:      vectors that are the same types
 69:   */
 70:   DMCreateGlobalVector(da,&x);
 71:   VecDuplicate(x,&r);

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

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

 91:   SNESSetFromOptions(snes);
 92:   SNESSolve(snes,NULL,x);

 94:   VecDestroy(&x);
 95:   VecDestroy(&r);
 96:   MatDestroy(&J);
 97:   MatNullSpaceDestroy(&constants);
 98:   SNESDestroy(&snes);
 99:   DMDestroy(&da);
100:   PetscFinalize();
101:   return ierr;
102: }

104: /*
105:    FormFunction - Evaluates nonlinear function, F(x).

107:    Input Parameters:
108: .  snes - the SNES context
109: .  x - input vector
110: .  ctx - optional user-defined context, as set by SNESSetFunction()

112:    Output Parameter:
113: .  f - function vector

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

129:   /* Get local work vector */
130:   DMGetLocalVector(da,&xlocal);

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

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

150:   /*
151:      Get local grid boundaries (for 1-dimensional DMDA):
152:        xs, xm  - starting grid index, width of local grid (no ghost points)
153:   */
154:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
155:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

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

164:   /*
165:      Restore vectors
166:   */
167:   DMDAVecRestoreArray(da,xlocal,&xx);
168:   DMDAVecRestoreArray(da,f,&ff);
169:   DMRestoreLocalVector(da,&xlocal);
170:   return(0);
171: }
172: /* ------------------------------------------------------------------- */
173: /*
174:    FormJacobian - Evaluates Jacobian matrix.

176:    Input Parameters:
177: .  snes - the SNES context
178: .  x - input vector
179: .  dummy - optional user-defined context (not used here)

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

196:   /*
197:      Get pointer to vector data
198:   */
199:   DMDAVecGetArrayRead(da,x,&xx);
200:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

202:   /*
203:     Get range of locally owned matrix
204:   */
205:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

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

219:   DMDAVecRestoreArrayRead(da,x,&xx);
220:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
221:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
222:   return(0);
223: }


226: /*TEST

228:    test:
229:       args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
230:       requires: !single

232: TEST*/