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 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 93 PetscScalar funcU(DMDACoor2d *coords) 94 { 95 return coords->x + coords->y; 96 } 97 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 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 */ 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 */ 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