1 static char help[] = "Surface processes in geophysics.\n\n";
2
3 #include <petscsnes.h>
4 #include <petscdm.h>
5 #include <petscdmda.h>
6
7 /*
8 User-defined application context - contains data needed by the
9 application-provided call-back routines, FormJacobianLocal() and
10 FormFunctionLocal().
11 */
12 typedef struct {
13 PetscReal D; /* The diffusion coefficient */
14 PetscReal K; /* The advection coefficient */
15 PetscInt m; /* Exponent for A */
16 } AppCtx;
17
18 /*
19 User-defined routines
20 */
21 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
22 extern PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscScalar **, Mat, AppCtx *);
23
main(int argc,char ** argv)24 int main(int argc, char **argv)
25 {
26 SNES snes; /* nonlinear solver */
27 AppCtx user; /* user-defined work context */
28 PetscInt its; /* iterations for convergence */
29 DM da;
30
31 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32 Initialize program
33 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37
38 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39 Initialize problem parameters
40 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Surface Process Problem Options", "SNES");
42 user.D = 1.0;
43 PetscCall(PetscOptionsReal("-D", "The diffusion coefficient D", __FILE__, user.D, &user.D, NULL));
44 user.K = 1.0;
45 PetscCall(PetscOptionsReal("-K", "The advection coefficient K", __FILE__, user.K, &user.K, NULL));
46 user.m = 1;
47 PetscCall(PetscOptionsInt("-m", "The exponent for A", __FILE__, user.m, &user.m, NULL));
48 PetscOptionsEnd();
49
50 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51 Create distributed array (DMDA) to manage parallel grid and vectors
52 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
54 PetscCall(DMSetFromOptions(da));
55 PetscCall(DMSetUp(da));
56 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
57 PetscCall(DMSetApplicationContext(da, &user));
58 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
59 PetscCall(SNESSetDM(snes, da));
60
61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62 Set local function evaluation routine
63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64 PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
65
66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67 Customize solver; set runtime options
68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69 PetscCall(SNESSetFromOptions(snes));
70
71 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 Solve nonlinear system
73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74 PetscCall(SNESSolve(snes, 0, 0));
75 PetscCall(SNESGetIterationNumber(snes, &its));
76
77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
80
81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 Free work space. All PETSc objects should be destroyed when they
83 are no longer needed.
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85
86 PetscCall(SNESDestroy(&snes));
87 PetscCall(DMDestroy(&da));
88
89 PetscCall(PetscFinalize());
90 return 0;
91 }
92
funcU(DMDACoor2d * coords)93 PetscScalar funcU(DMDACoor2d *coords)
94 {
95 return coords->x + coords->y;
96 }
97
funcA(PetscScalar z,AppCtx * user)98 PetscScalar funcA(PetscScalar z, AppCtx *user)
99 {
100 PetscScalar v = 1.0;
101 PetscInt i;
102
103 for (i = 0; i < user->m; ++i) v *= z;
104 return v;
105 }
106
funcADer(PetscScalar z,AppCtx * user)107 PetscScalar funcADer(PetscScalar z, AppCtx *user)
108 {
109 PetscScalar v = 1.0;
110 PetscInt i;
111
112 for (i = 0; i < user->m - 1; ++i) v *= z;
113 return (PetscScalar)user->m * v;
114 }
115
116 /*
117 FormFunctionLocal - Evaluates nonlinear function, F(x).
118 */
FormFunctionLocal(DMDALocalInfo * info,PetscScalar ** x,PetscScalar ** f,AppCtx * user)119 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
120 {
121 DM coordDA;
122 Vec coordinates;
123 DMDACoor2d **coords;
124 PetscScalar u, ux, uy, uxx, uyy;
125 PetscReal D, K, hx, hy, hxdhy, hydhx;
126 PetscInt i, j;
127
128 PetscFunctionBeginUser;
129 D = user->D;
130 K = user->K;
131 hx = 1.0 / (PetscReal)(info->mx - 1);
132 hy = 1.0 / (PetscReal)(info->my - 1);
133 hxdhy = hx / hy;
134 hydhx = hy / hx;
135 /*
136 Compute function over the locally owned part of the grid
137 */
138 PetscCall(DMGetCoordinateDM(info->da, &coordDA));
139 PetscCall(DMGetCoordinates(info->da, &coordinates));
140 PetscCall(DMDAVecGetArray(coordDA, coordinates, &coords));
141 for (j = info->ys; j < info->ys + info->ym; j++) {
142 for (i = info->xs; i < info->xs + info->xm; i++) {
143 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) f[j][i] = x[j][i];
144 else {
145 u = x[j][i];
146 ux = (x[j][i + 1] - x[j][i]) / hx;
147 uy = (x[j + 1][i] - x[j][i]) / hy;
148 uxx = (2.0 * u - x[j][i - 1] - x[j][i + 1]) * hydhx;
149 uyy = (2.0 * u - x[j - 1][i] - x[j + 1][i]) * hxdhy;
150 f[j][i] = D * (uxx + uyy) - (K * funcA(x[j][i], user) * PetscSqrtScalar(ux * ux + uy * uy) + funcU(&coords[j][i])) * hx * hy;
151 PetscCheck(!PetscIsInfOrNanScalar(f[j][i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Invalid residual: %g", (double)PetscRealPart(f[j][i]));
152 }
153 }
154 }
155 PetscCall(DMDAVecRestoreArray(coordDA, coordinates, &coords));
156 PetscCall(PetscLogFlops(11.0 * info->ym * info->xm));
157 PetscFunctionReturn(PETSC_SUCCESS);
158 }
159
160 /*
161 FormJacobianLocal - Evaluates Jacobian matrix.
162 */
FormJacobianLocal(DMDALocalInfo * info,PetscScalar ** x,Mat jac,AppCtx * user)163 PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscScalar **x, Mat jac, AppCtx *user)
164 {
165 MatStencil col[5], row;
166 PetscScalar D, K, A, v[5], hx, hy, hxdhy, hydhx, ux, uy;
167 PetscReal normGradZ;
168 PetscInt i, j, k;
169
170 PetscFunctionBeginUser;
171 D = user->D;
172 K = user->K;
173 hx = 1.0 / (PetscReal)(info->mx - 1);
174 hy = 1.0 / (PetscReal)(info->my - 1);
175 hxdhy = hx / hy;
176 hydhx = hy / hx;
177
178 /*
179 Compute entries for the locally owned part of the Jacobian.
180 - Currently, all PETSc parallel matrix formats are partitioned by
181 contiguous chunks of rows across the processors.
182 - Each processor needs to insert only elements that it owns
183 locally (but any non-local elements will be sent to the
184 appropriate processor during matrix assembly).
185 - Here, we set all entries for a particular row at once.
186 - We can set matrix entries either using either
187 MatSetValuesLocal() or MatSetValues(), as discussed above.
188 */
189 for (j = info->ys; j < info->ys + info->ym; j++) {
190 for (i = info->xs; i < info->xs + info->xm; i++) {
191 row.j = j;
192 row.i = i;
193 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
194 /* boundary points */
195 v[0] = 1.0;
196 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
197 } else {
198 /* interior grid points */
199 ux = (x[j][i + 1] - x[j][i]) / hx;
200 uy = (x[j + 1][i] - x[j][i]) / hy;
201 normGradZ = PetscRealPart(PetscSqrtScalar(ux * ux + uy * uy));
202 if (normGradZ < 1.0e-8) normGradZ = 1.0e-8;
203 A = funcA(x[j][i], user);
204
205 v[0] = -D * hxdhy;
206 col[0].j = j - 1;
207 col[0].i = i;
208 v[1] = -D * hydhx;
209 col[1].j = j;
210 col[1].i = i - 1;
211 v[2] = D * 2.0 * (hydhx + hxdhy) + K * (funcADer(x[j][i], user) * normGradZ - A / normGradZ) * hx * hy;
212 col[2].j = row.j;
213 col[2].i = row.i;
214 v[3] = -D * hydhx + K * A * hx * hy / (2.0 * normGradZ);
215 col[3].j = j;
216 col[3].i = i + 1;
217 v[4] = -D * hxdhy + K * A * hx * hy / (2.0 * normGradZ);
218 col[4].j = j + 1;
219 col[4].i = i;
220 for (k = 0; k < 5; ++k) PetscCheck(!PetscIsInfOrNanScalar(v[k]), PETSC_COMM_SELF, PETSC_ERR_FP, "Invalid residual: %g", (double)PetscRealPart(v[k]));
221 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
222 }
223 }
224 }
225
226 /*
227 Assemble matrix, using the 2-step process:
228 MatAssemblyBegin(), MatAssemblyEnd().
229 */
230 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
231 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
232 /*
233 Tell the matrix we will never add a new nonzero location to the
234 matrix. If we do, it will generate an error.
235 */
236 PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
237 PetscFunctionReturn(PETSC_SUCCESS);
238 }
239
240 /*TEST
241
242 test:
243 args: -snes_view -snes_monitor_short -da_refine 1 -pc_type mg -ksp_type fgmres -pc_mg_type full -mg_levels_ksp_chebyshev_esteig 0.5,1.1
244
245 test:
246 suffix: ew_1
247 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 1
248 requires: !single
249
250 test:
251 suffix: ew_2
252 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 2
253
254 test:
255 suffix: ew_3
256 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 3
257 requires: !single
258
259 test:
260 suffix: fm_rise_2
261 args: -K 3 -m 1 -D 0.2 -snes_monitor_short -snes_type ngmres -snes_npc_side right -npc_snes_type newtonls -npc_snes_linesearch_type basic -snes_ngmres_restart_it 1 -snes_ngmres_restart_fm_rise
262 requires: !single
263
264 test:
265 suffix: fm_rise_4
266 args: -K 3 -m 1 -D 0.2 -snes_monitor_short -snes_type ngmres -snes_npc_side right -npc_snes_type newtonls -npc_snes_linesearch_type basic -snes_ngmres_restart_it 2 -snes_ngmres_restart_fm_rise -snes_rtol 1.e-2 -snes_max_it 5
267
268 TEST*/
269