1 static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n"; 2 3 /* 4 The linear and nonlinear versions of these should give almost identical results on this problem 5 6 Richardson 7 Nonlinear: 8 -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor 9 10 Linear: 11 -snes_rtol 1.e-12 -snes_monitor -ksp_rtol 1.e-12 -ksp_monitor -ksp_type richardson -pc_type none -ksp_richardson_self_scale -info 12 13 GMRES 14 Nonlinear: 15 -snes_rtol 1.e-12 -snes_monitor -snes_type ngmres 16 17 Linear: 18 -snes_rtol 1.e-12 -snes_monitor -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none 19 20 CG 21 Nonlinear: 22 -snes_rtol 1.e-12 -snes_monitor -snes_type ncg -snes_linesearch_monitor 23 24 Linear: 25 -snes_rtol 1.e-12 -snes_monitor -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none 26 27 Multigrid 28 Linear: 29 1 level: 30 -snes_rtol 1.e-12 -snes_monitor -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor 31 -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12 -ksp_monitor_true_residual 32 33 n levels: 34 -da_refine n 35 36 Nonlinear: 37 1 level: 38 -snes_rtol 1.e-12 -snes_monitor -snes_type fas -fas_levels_snes_monitor 39 40 n levels: 41 -da_refine n -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly 42 */ 43 44 /* 45 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 46 Include "petscsnes.h" so that we can use SNES solvers. Note that this 47 */ 48 #include <petscdm.h> 49 #include <petscdmda.h> 50 #include <petscsnes.h> 51 52 /* 53 User-defined routines 54 */ 55 extern PetscErrorCode FormMatrix(DM, Mat); 56 extern PetscErrorCode MyComputeFunction(SNES, Vec, Vec, void *); 57 extern PetscErrorCode MyComputeJacobian(SNES, Vec, Mat, Mat, void *); 58 extern PetscErrorCode NonlinearGS(SNES, Vec); 59 60 int main(int argc, char **argv) 61 { 62 SNES snes; /* nonlinear solver */ 63 SNES psnes; /* nonlinear Gauss-Seidel approximate solver */ 64 Vec x, b; /* solution vector */ 65 PetscInt its; /* iterations for convergence */ 66 DM da; 67 PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */ 68 69 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 70 Initialize program 71 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 72 73 PetscFunctionBeginUser; 74 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 75 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Create nonlinear solver context 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 79 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 80 81 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_ngs_as_npc", &use_ngs_as_npc, 0)); 82 83 if (use_ngs_as_npc) { 84 PetscCall(SNESGetNPC(snes, &psnes)); 85 PetscCall(SNESSetType(psnes, SNESSHELL)); 86 PetscCall(SNESShellSetSolve(psnes, NonlinearGS)); 87 } 88 89 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 90 Create distributed array (DMDA) to manage parallel grid and vectors 91 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 92 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)); 93 PetscCall(DMSetFromOptions(da)); 94 PetscCall(DMSetUp(da)); 95 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 96 PetscCall(SNESSetDM(snes, da)); 97 if (use_ngs_as_npc) PetscCall(SNESShellSetContext(psnes, da)); 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Extract global vectors from DMDA; then duplicate for remaining 100 vectors that are the same types 101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 102 PetscCall(DMCreateGlobalVector(da, &x)); 103 PetscCall(DMCreateGlobalVector(da, &b)); 104 PetscCall(VecSet(b, 1.0)); 105 106 PetscCall(SNESSetFunction(snes, NULL, MyComputeFunction, NULL)); 107 PetscCall(SNESSetJacobian(snes, NULL, NULL, MyComputeJacobian, NULL)); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Customize nonlinear solver; set runtime options 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 PetscCall(SNESSetFromOptions(snes)); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Solve nonlinear system 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 PetscCall(SNESSolve(snes, b, x)); 118 PetscCall(SNESGetIterationNumber(snes, &its)); 119 120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 Free work space. All PETSc objects should be destroyed when they 125 are no longer needed. 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 PetscCall(VecDestroy(&x)); 128 PetscCall(VecDestroy(&b)); 129 PetscCall(SNESDestroy(&snes)); 130 PetscCall(DMDestroy(&da)); 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 /* ------------------------------------------------------------------- */ 136 PetscErrorCode MyComputeFunction(SNES snes, Vec x, Vec F, PetscCtx ctx) 137 { 138 Mat J; 139 DM dm; 140 141 PetscFunctionBeginUser; 142 PetscCall(SNESGetDM(snes, &dm)); 143 PetscCall(PetscObjectQuery((PetscObject)dm, "_ex35_J", (PetscObject *)&J)); 144 if (!J) { 145 PetscCall(DMSetMatType(dm, MATAIJ)); 146 PetscCall(DMCreateMatrix(dm, &J)); 147 PetscCall(MatSetDM(J, NULL)); 148 PetscCall(FormMatrix(dm, J)); 149 PetscCall(PetscObjectCompose((PetscObject)dm, "_ex35_J", (PetscObject)J)); 150 PetscCall(PetscObjectDereference((PetscObject)J)); 151 } 152 PetscCall(MatMult(J, x, F)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 PetscErrorCode MyComputeJacobian(SNES snes, Vec x, Mat J, Mat Jp, PetscCtx ctx) 157 { 158 DM dm; 159 160 PetscFunctionBeginUser; 161 PetscCall(SNESGetDM(snes, &dm)); 162 PetscCall(FormMatrix(dm, Jp)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 PetscErrorCode FormMatrix(DM da, Mat jac) 167 { 168 PetscInt i, j, nrows = 0; 169 MatStencil col[5], row, *rows; 170 PetscScalar v[5], hx, hy, hxdhy, hydhx; 171 DMDALocalInfo info; 172 173 PetscFunctionBeginUser; 174 PetscCall(DMDAGetLocalInfo(da, &info)); 175 hx = 1.0 / (PetscReal)(info.mx - 1); 176 hy = 1.0 / (PetscReal)(info.my - 1); 177 hxdhy = hx / hy; 178 hydhx = hy / hx; 179 180 PetscCall(PetscMalloc1(info.ym * info.xm, &rows)); 181 /* 182 Compute entries for the locally owned part of the Jacobian. 183 - Currently, all PETSc parallel matrix formats are partitioned by 184 contiguous chunks of rows across the processors. 185 - Each processor needs to insert only elements that it owns 186 locally (but any non-local elements will be sent to the 187 appropriate processor during matrix assembly). 188 - Here, we set all entries for a particular row at once. 189 - We can set matrix entries either using either 190 MatSetValuesLocal() or MatSetValues(), as discussed above. 191 */ 192 for (j = info.ys; j < info.ys + info.ym; j++) { 193 for (i = info.xs; i < info.xs + info.xm; i++) { 194 row.j = j; 195 row.i = i; 196 /* boundary points */ 197 if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) { 198 v[0] = 2.0 * (hydhx + hxdhy); 199 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES)); 200 rows[nrows].i = i; 201 rows[nrows++].j = j; 202 } else { 203 /* interior grid points */ 204 v[0] = -hxdhy; 205 col[0].j = j - 1; 206 col[0].i = i; 207 v[1] = -hydhx; 208 col[1].j = j; 209 col[1].i = i - 1; 210 v[2] = 2.0 * (hydhx + hxdhy); 211 col[2].j = row.j; 212 col[2].i = row.i; 213 v[3] = -hydhx; 214 col[3].j = j; 215 col[3].i = i + 1; 216 v[4] = -hxdhy; 217 col[4].j = j + 1; 218 col[4].i = i; 219 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES)); 220 } 221 } 222 } 223 224 /* 225 Assemble matrix, using the 2-step process: 226 MatAssemblyBegin(), MatAssemblyEnd(). 227 */ 228 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 229 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 230 PetscCall(MatZeroRowsColumnsStencil(jac, nrows, rows, 2.0 * (hydhx + hxdhy), NULL, NULL)); 231 PetscCall(PetscFree(rows)); 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 /* ------------------------------------------------------------------- */ 241 /* 242 Applies some sweeps on nonlinear Gauss-Seidel on each process 243 244 */ 245 PetscErrorCode NonlinearGS(SNES snes, Vec X) 246 { 247 PetscInt i, j, Mx, My, xs, ys, xm, ym, its, l; 248 PetscReal hx, hy, hxdhy, hydhx; 249 PetscScalar **x, F, J, u, uxx, uyy; 250 DM da; 251 Vec localX; 252 253 PetscFunctionBeginUser; 254 PetscCall(SNESGetTolerances(snes, NULL, NULL, NULL, &its, NULL)); 255 PetscCall(SNESShellGetContext(snes, &da)); 256 257 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 258 259 hx = 1.0 / (PetscReal)(Mx - 1); 260 hy = 1.0 / (PetscReal)(My - 1); 261 hxdhy = hx / hy; 262 hydhx = hy / hx; 263 264 PetscCall(DMGetLocalVector(da, &localX)); 265 266 for (l = 0; l < its; l++) { 267 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 268 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 269 /* 270 Get a pointer to vector data. 271 - For default PETSc vectors, VecGetArray() returns a pointer to 272 the data array. Otherwise, the routine is implementation dependent. 273 - You MUST call VecRestoreArray() when you no longer need access to 274 the array. 275 */ 276 PetscCall(DMDAVecGetArray(da, localX, &x)); 277 278 /* 279 Get local grid boundaries (for 2-dimensional DMDA): 280 xs, ys - starting grid indices (no ghost points) 281 xm, ym - widths of local grid (no ghost points) 282 283 */ 284 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 285 286 for (j = ys; j < ys + ym; j++) { 287 for (i = xs; i < xs + xm; i++) { 288 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) { 289 /* boundary conditions are all zero Dirichlet */ 290 x[j][i] = 0.0; 291 } else { 292 u = x[j][i]; 293 uxx = (2.0 * u - x[j][i - 1] - x[j][i + 1]) * hydhx; 294 uyy = (2.0 * u - x[j - 1][i] - x[j + 1][i]) * hxdhy; 295 F = uxx + uyy; 296 J = 2.0 * (hydhx + hxdhy); 297 u = u - F / J; 298 299 x[j][i] = u; 300 } 301 } 302 } 303 304 /* 305 Restore vector 306 */ 307 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 308 PetscCall(DMLocalToGlobalBegin(da, localX, INSERT_VALUES, X)); 309 PetscCall(DMLocalToGlobalEnd(da, localX, INSERT_VALUES, X)); 310 } 311 PetscCall(DMRestoreLocalVector(da, &localX)); 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 /*TEST 316 317 test: 318 args: -snes_monitor_short -snes_type nrichardson 319 requires: !single 320 321 test: 322 suffix: 2 323 args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale 324 requires: !single 325 326 test: 327 suffix: 3 328 args: -snes_monitor_short -snes_type ngmres 329 330 test: 331 suffix: 4 332 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none 333 334 test: 335 suffix: 5 336 args: -snes_monitor_short -snes_type ncg 337 338 test: 339 suffix: 6 340 args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none 341 342 test: 343 suffix: 7 344 args: -da_refine 2 -snes_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor_short -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor_short 345 requires: !single 346 347 test: 348 suffix: 8 349 args: -da_refine 2 -snes_monitor_short -snes_type fas -fas_levels_snes_monitor_short -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_type fas -snes_rtol 1.e-5 350 351 test: 352 suffix: 9 353 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc 354 355 test: 356 suffix: 10 357 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false 358 359 TEST*/ 360