Actual source code: ex46.c
1: static char help[] = "Surface processes in geophysics.\n\n";
3: #include <petscsnes.h>
4: #include <petscdm.h>
5: #include <petscdmda.h>
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;
18: /*
19: User-defined routines
20: */
21: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
22: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo *, PetscScalar **, Mat, AppCtx *);
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;
31: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: Initialize program
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
35: PetscFunctionBeginUser;
36: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
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();
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));
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Set local function evaluation routine
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Customize solver; set runtime options
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: PetscCall(SNESSetFromOptions(snes));
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Solve nonlinear system
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscCall(SNESSolve(snes, 0, 0));
75: PetscCall(SNESGetIterationNumber(snes, &its));
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Free work space. All PETSc objects should be destroyed when they
83: are no longer needed.
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: PetscCall(SNESDestroy(&snes));
87: PetscCall(DMDestroy(&da));
89: PetscCall(PetscFinalize());
90: return 0;
91: }
93: PetscScalar funcU(DMDACoor2d *coords)
94: {
95: return coords->x + coords->y;
96: }
98: PetscScalar funcA(PetscScalar z, AppCtx *user)
99: {
100: PetscScalar v = 1.0;
101: PetscInt i;
103: for (i = 0; i < user->m; ++i) v *= z;
104: return v;
105: }
107: PetscScalar funcADer(PetscScalar z, AppCtx *user)
108: {
109: PetscScalar v = 1.0;
110: PetscInt i;
112: for (i = 0; i < user->m - 1; ++i) v *= z;
113: return (PetscScalar)user->m * v;
114: }
116: /*
117: FormFunctionLocal - Evaluates nonlinear function, F(x).
118: */
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;
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: }
160: /*
161: FormJacobianLocal - Evaluates Jacobian matrix.
162: */
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;
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;
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);
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: }
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: }
240: /*TEST
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
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
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
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
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
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
268: TEST*/