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