1 /* ------------------------------------------------------------------------
2
3 Solid Fuel Ignition (SFI) problem. This problem is modeled by the
4 partial differential equation
5
6 -Laplacian(u) - lambda * exp(u) = 0, 0 < x,y,z < 1,
7
8 with boundary conditions
9
10 u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
11
12 A finite difference approximation with the usual 7-point stencil
13 is used to discretize the boundary value problem to obtain a
14 nonlinear system of equations. The problem is solved in a 3D
15 rectangular domain, using distributed arrays (DAs) to partition
16 the parallel grid.
17
18 ------------------------------------------------------------------------- */
19
20 #include "Bratu3Dimpl.h"
21
FormInitGuess(DM da,Vec X,Params * p)22 PetscErrorCode FormInitGuess(DM da, Vec X, Params *p)
23 {
24 PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
25 PetscReal lambda,temp1,hx,hy,hz,tempk,tempj;
26 PetscScalar ***x;
27
28 PetscFunctionBegin;
29 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
30
31 lambda = p->lambda_;
32 hx = 1.0/(PetscReal)(Mx-1);
33 hy = 1.0/(PetscReal)(My-1);
34 hz = 1.0/(PetscReal)(Mz-1);
35 temp1 = lambda/(lambda + 1.0);
36
37 /*
38 Get a pointer to vector data.
39
40 - For default PETSc vectors, VecGetArray() returns a pointer to
41 the data array. Otherwise, the routine is implementation
42 dependent.
43
44 - You MUST call VecRestoreArray() when you no longer need access
45 to the array.
46 */
47 PetscCall(DMDAVecGetArray(da,X,&x));
48
49 /*
50 Get local grid boundaries (for 3-dimensional DMDA):
51
52 - xs, ys, zs: starting grid indices (no ghost points)
53
54 - xm, ym, zm: widths of local grid (no ghost points)
55 */
56 PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
57
58 /*
59 Compute initial guess over the locally owned part of the grid
60 */
61 for (k=zs; k<zs+zm; k++) {
62 tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
63 for (j=ys; j<ys+ym; j++) {
64 tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
65 for (i=xs; i<xs+xm; i++) {
66 if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
67 /* boundary conditions are all zero Dirichlet */
68 x[k][j][i] = 0.0;
69 } else {
70 x[k][j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
71 }
72 }
73 }
74 }
75
76 /*
77 Restore vector
78 */
79 PetscCall(DMDAVecRestoreArray(da,X,&x));
80 PetscFunctionReturn(PETSC_SUCCESS);
81 }
82
FormFunction(DM da,Vec X,Vec F,Params * p)83 PetscErrorCode FormFunction(DM da, Vec X, Vec F, Params *p)
84 {
85 PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
86 PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
87 PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
88 Vec localX;
89
90 PetscFunctionBegin;
91 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
92 lambda = p->lambda_;
93 hx = 1.0/(PetscReal)(Mx-1);
94 hy = 1.0/(PetscReal)(My-1);
95 hz = 1.0/(PetscReal)(Mz-1);
96 sc = hx*hy*hz*lambda;
97 hxhzdhy = hx*hz/hy;
98 hyhzdhx = hy*hz/hx;
99 hxhydhz = hx*hy/hz;
100
101 /*
102
103 */
104 PetscCall(DMGetLocalVector(da,&localX));
105
106 /*
107 Scatter ghost points to local vector,using the 2-step process
108 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). By placing code
109 between these two statements, computations can be done while
110 messages are in transition.
111 */
112 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX));
113 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX));
114
115 /*
116 Get pointers to vector data.
117 */
118 PetscCall(DMDAVecGetArray(da,localX,&x));
119 PetscCall(DMDAVecGetArray(da,F,&f));
120
121 /*
122 Get local grid boundaries.
123 */
124 PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
125
126 /*
127 Compute function over the locally owned part of the grid.
128 */
129 for (k=zs; k<zs+zm; k++) {
130 for (j=ys; j<ys+ym; j++) {
131 for (i=xs; i<xs+xm; i++) {
132 if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
133 /* boundary points */
134 f[k][j][i] = x[k][j][i] - 0.0;
135 } else {
136 /* interior grid points */
137 u = x[k][j][i];
138 u_east = x[k][j][i+1];
139 u_west = x[k][j][i-1];
140 u_north = x[k][j+1][i];
141 u_south = x[k][j-1][i];
142 u_up = x[k+1][j][i];
143 u_down = x[k-1][j][i];
144 u_xx = (-u_east + two*u - u_west)*hyhzdhx;
145 u_yy = (-u_north + two*u - u_south)*hxhzdhy;
146 u_zz = (-u_up + two*u - u_down)*hxhydhz;
147 f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
148 }
149 }
150 }
151 }
152
153 /*
154 Restore vectors.
155 */
156 PetscCall(DMDAVecRestoreArray(da,F,&f));
157 PetscCall(DMDAVecRestoreArray(da,localX,&x));
158 PetscCall(DMRestoreLocalVector(da,&localX));
159 PetscCall(PetscLogFlops(11.0*ym*xm));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
FormJacobian(DM da,Vec X,Mat J,Params * p)163 PetscErrorCode FormJacobian(DM da, Vec X, Mat J, Params *p)
164 {
165 PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
166 PetscReal lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
167 PetscScalar v[7],***x;
168 MatStencil col[7],row;
169 Vec localX;
170
171 PetscFunctionBegin;
172 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE));
173 lambda = p->lambda_;
174 hx = 1.0/(PetscReal)(Mx-1);
175 hy = 1.0/(PetscReal)(My-1);
176 hz = 1.0/(PetscReal)(Mz-1);
177 sc = hx*hy*hz*lambda;
178 hxhzdhy = hx*hz/hy;
179 hyhzdhx = hy*hz/hx;
180 hxhydhz = hx*hy/hz;
181
182 /*
183
184 */
185 PetscCall(DMGetLocalVector(da,&localX));
186
187 /*
188 Scatter ghost points to local vector, using the 2-step process
189 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). By placing code
190 between these two statements, computations can be done while
191 messages are in transition.
192 */
193 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX));
194 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX));
195
196 /*
197 Get pointer to vector data.
198 */
199 PetscCall(DMDAVecGetArray(da,localX,&x));
200
201 /*
202 Get local grid boundaries.
203 */
204 PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
205
206 /*
207 Compute entries for the locally owned part of the Jacobian.
208
209 - Currently, all PETSc parallel matrix formats are partitioned by
210 contiguous chunks of rows across the processors.
211
212 - Each processor needs to insert only elements that it owns
213 locally (but any non-local elements will be sent to the
214 appropriate processor during matrix assembly).
215
216 - Here, we set all entries for a particular row at once.
217
218 - We can set matrix entries either using either
219 MatSetValuesLocal() or MatSetValues(), as discussed above.
220 */
221 for (k=zs; k<zs+zm; k++) {
222 for (j=ys; j<ys+ym; j++) {
223 for (i=xs; i<xs+xm; i++) {
224 row.k = k; row.j = j; row.i = i;
225 /* boundary points */
226 if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
227 v[0] = 1.0;
228 PetscCall(MatSetValuesStencil(J,1,&row,1,&row,v,INSERT_VALUES));
229 } else {
230 /* interior grid points */
231 v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i;
232 v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i;
233 v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1;
234 v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
235 v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1;
236 v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i;
237 v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i;
238 PetscCall(MatSetValuesStencil(J,1,&row,7,col,v,INSERT_VALUES));
239 }
240 }
241 }
242 }
243 PetscCall(DMDAVecRestoreArray(da,localX,&x));
244 PetscCall(DMRestoreLocalVector(da,&localX));
245
246 /*
247 Assemble matrix, using the 2-step process: MatAssemblyBegin(),
248 MatAssemblyEnd().
249 */
250 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY));
251 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY));
252
253 /*
254 Tell the matrix we will never add a new nonzero location to the
255 matrix. If we do, it will generate an error.
256 */
257 PetscCall(MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE));
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260