1
2 static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n";
3
4 /*T
5 Concepts: SNES^basic parallel example
6 Concepts: periodic boundary conditions
7 Processors: n
8 T*/
9
10
11
12 /*
13 Compare this example to ex3.c that handles Dirichlet boundary conditions
14
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.
17
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 */
27
28 #include <petscdm.h>
29 #include <petscdmda.h>
30 #include <petscsnes.h>
31
32
33 PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
34 PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
35
main(int argc,char ** argv)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 */
42 PetscErrorCode ierr;
43 PetscInt N = 5;
44 MatNullSpace constants;
45
46 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
47 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);CHKERRQ(ierr);
48
49 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50 Create nonlinear solver context
51 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52
53 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
54
55 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56 Create vector data structures; set function evaluation routine
57 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58
59 /*
60 Create distributed array (DMDA) to manage parallel grid and vectors
61 */
62 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);CHKERRQ(ierr);
63 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
64 ierr = DMSetUp(da);CHKERRQ(ierr);
65
66 /*
67 Extract global and local vectors from DMDA; then duplicate for remaining
68 vectors that are the same types
69 */
70 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
71 ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
72
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 ierr = SNESSetFunction(snes,r,FormFunction,da);CHKERRQ(ierr);
82
83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84 Create matrix data structure; set Jacobian evaluation routine
85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
87 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);CHKERRQ(ierr);
88 ierr = MatSetNullSpace(J,constants);CHKERRQ(ierr);
89 ierr = SNESSetJacobian(snes,J,J,FormJacobian,da);CHKERRQ(ierr);
90
91 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
92 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
93
94 ierr = VecDestroy(&x);CHKERRQ(ierr);
95 ierr = VecDestroy(&r);CHKERRQ(ierr);
96 ierr = MatDestroy(&J);CHKERRQ(ierr);
97 ierr = MatNullSpaceDestroy(&constants);CHKERRQ(ierr);
98 ierr = SNESDestroy(&snes);CHKERRQ(ierr);
99 ierr = DMDestroy(&da);CHKERRQ(ierr);
100 ierr = PetscFinalize();
101 return ierr;
102 }
103
104 /*
105 FormFunction - Evaluates nonlinear function, F(x).
106
107 Input Parameters:
108 . snes - the SNES context
109 . x - input vector
110 . ctx - optional user-defined context, as set by SNESSetFunction()
111
112 Output Parameter:
113 . f - function vector
114
115 Note:
116 The user-defined context can contain any application-specific
117 data needed for the function evaluation.
118 */
FormFunction(SNES snes,Vec x,Vec f,void * ctx)119 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
120 {
121 DM da = (DM) ctx;
122 PetscScalar *xx,*ff;
123 PetscReal h;
124 PetscErrorCode ierr;
125 PetscInt i,M,xs,xm;
126 Vec xlocal;
127
128 PetscFunctionBeginUser;
129 /* Get local work vector */
130 ierr = DMGetLocalVector(da,&xlocal);CHKERRQ(ierr);
131
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 ierr = DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);
139 ierr = DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);CHKERRQ(ierr);
140
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 ierr = DMDAVecGetArray(da,xlocal,&xx);CHKERRQ(ierr);
148 ierr = DMDAVecGetArray(da,f,&ff);CHKERRQ(ierr);
149
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 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
155 ierr = DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
156
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);
163
164 /*
165 Restore vectors
166 */
167 ierr = DMDAVecRestoreArray(da,xlocal,&xx);CHKERRQ(ierr);
168 ierr = DMDAVecRestoreArray(da,f,&ff);CHKERRQ(ierr);
169 ierr = DMRestoreLocalVector(da,&xlocal);CHKERRQ(ierr);
170 PetscFunctionReturn(0);
171 }
172 /* ------------------------------------------------------------------- */
173 /*
174 FormJacobian - Evaluates Jacobian matrix.
175
176 Input Parameters:
177 . snes - the SNES context
178 . x - input vector
179 . dummy - optional user-defined context (not used here)
180
181 Output Parameters:
182 . jac - Jacobian matrix
183 . B - optionally different preconditioning matrix
184 . flag - flag indicating matrix structure
185 */
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * ctx)186 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
187 {
188 PetscScalar *xx,A[3];
189 PetscErrorCode ierr;
190 PetscInt i,M,xs,xm;
191 DM da = (DM) ctx;
192 MatStencil row,cols[3];
193 PetscReal h;
194
195 PetscFunctionBeginUser;
196 /*
197 Get pointer to vector data
198 */
199 ierr = DMDAVecGetArrayRead(da,x,&xx);CHKERRQ(ierr);
200 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
201
202 /*
203 Get range of locally owned matrix
204 */
205 ierr = DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
206
207 ierr = MatZeroEntries(jac);CHKERRQ(ierr);
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 ierr = MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);CHKERRQ(ierr);
217 }
218
219 ierr = DMDAVecRestoreArrayRead(da,x,&xx);CHKERRQ(ierr);
220 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
221 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
222 PetscFunctionReturn(0);
223 }
224
225
226 /*TEST
227
228 test:
229 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
230 requires: !single
231
232 TEST*/
233