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:   */
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);
217:   }

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