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 { 64 SNES snes; /* nonlinear solver */ 65 SNES psnes; /* nonlinear Gauss-Seidel approximate solver */ 66 Vec x,b; /* solution vector */ 67 PetscInt its; /* iterations for convergence */ 68 DM da; 69 PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */ 70 71 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 72 Initialize program 73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 74 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) { 99 PetscCall(SNESShellSetContext(psnes,da)); 100 } 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Extract global vectors from DMDA; then duplicate for remaining 103 vectors that are the same types 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 PetscCall(DMCreateGlobalVector(da,&x)); 106 PetscCall(DMCreateGlobalVector(da,&b)); 107 PetscCall(VecSet(b,1.0)); 108 109 PetscCall(SNESSetFunction(snes,NULL,MyComputeFunction,NULL)); 110 PetscCall(SNESSetJacobian(snes,NULL,NULL,MyComputeJacobian,NULL)); 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Customize nonlinear solver; set runtime options 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 PetscCall(SNESSetFromOptions(snes)); 116 117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118 Solve nonlinear system 119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 120 PetscCall(SNESSolve(snes,b,x)); 121 PetscCall(SNESGetIterationNumber(snes,&its)); 122 123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Free work space. All PETSc objects should be destroyed when they 128 are no longer needed. 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 PetscCall(VecDestroy(&x)); 131 PetscCall(VecDestroy(&b)); 132 PetscCall(SNESDestroy(&snes)); 133 PetscCall(DMDestroy(&da)); 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 138 /* ------------------------------------------------------------------- */ 139 PetscErrorCode MyComputeFunction(SNES snes,Vec x,Vec F,void *ctx) 140 { 141 Mat J; 142 DM dm; 143 144 PetscFunctionBeginUser; 145 PetscCall(SNESGetDM(snes,&dm)); 146 PetscCall(DMGetApplicationContext(dm,&J)); 147 if (!J) { 148 PetscCall(DMSetMatType(dm,MATAIJ)); 149 PetscCall(DMCreateMatrix(dm,&J)); 150 PetscCall(MatSetDM(J, NULL)); 151 PetscCall(FormMatrix(dm,J)); 152 PetscCall(DMSetApplicationContext(dm,J)); 153 PetscCall(DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy)); 154 } 155 PetscCall(MatMult(J,x,F)); 156 PetscFunctionReturn(0); 157 } 158 159 PetscErrorCode MyComputeJacobian(SNES snes,Vec x,Mat J,Mat Jp,void *ctx) 160 { 161 DM dm; 162 163 PetscFunctionBeginUser; 164 PetscCall(SNESGetDM(snes,&dm)); 165 PetscCall(FormMatrix(dm,Jp)); 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode FormMatrix(DM da,Mat jac) 170 { 171 PetscInt i,j,nrows = 0; 172 MatStencil col[5],row,*rows; 173 PetscScalar v[5],hx,hy,hxdhy,hydhx; 174 DMDALocalInfo info; 175 176 PetscFunctionBeginUser; 177 PetscCall(DMDAGetLocalInfo(da,&info)); 178 hx = 1.0/(PetscReal)(info.mx-1); 179 hy = 1.0/(PetscReal)(info.my-1); 180 hxdhy = hx/hy; 181 hydhx = hy/hx; 182 183 PetscCall(PetscMalloc1(info.ym*info.xm,&rows)); 184 /* 185 Compute entries for the locally owned part of the Jacobian. 186 - Currently, all PETSc parallel matrix formats are partitioned by 187 contiguous chunks of rows across the processors. 188 - Each processor needs to insert only elements that it owns 189 locally (but any non-local elements will be sent to the 190 appropriate processor during matrix assembly). 191 - Here, we set all entries for a particular row at once. 192 - We can set matrix entries either using either 193 MatSetValuesLocal() or MatSetValues(), as discussed above. 194 */ 195 for (j=info.ys; j<info.ys+info.ym; j++) { 196 for (i=info.xs; i<info.xs+info.xm; i++) { 197 row.j = j; row.i = i; 198 /* boundary points */ 199 if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) { 200 v[0] = 2.0*(hydhx + hxdhy); 201 PetscCall(MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES)); 202 rows[nrows].i = i; 203 rows[nrows++].j = j; 204 } else { 205 /* interior grid points */ 206 v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i; 207 v[1] = -hydhx; col[1].j = j; col[1].i = i-1; 208 v[2] = 2.0*(hydhx + hxdhy); col[2].j = row.j; col[2].i = row.i; 209 v[3] = -hydhx; col[3].j = j; col[3].i = i+1; 210 v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i; 211 PetscCall(MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES)); 212 } 213 } 214 } 215 216 /* 217 Assemble matrix, using the 2-step process: 218 MatAssemblyBegin(), MatAssemblyEnd(). 219 */ 220 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 221 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 222 PetscCall(MatZeroRowsColumnsStencil(jac,nrows,rows,2.0*(hydhx + hxdhy),NULL,NULL)); 223 PetscCall(PetscFree(rows)); 224 /* 225 Tell the matrix we will never add a new nonzero location to the 226 matrix. If we do, it will generate an error. 227 */ 228 PetscCall(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 229 PetscFunctionReturn(0); 230 } 231 232 /* ------------------------------------------------------------------- */ 233 /* 234 Applies some sweeps on nonlinear Gauss-Seidel on each process 235 236 */ 237 PetscErrorCode NonlinearGS(SNES snes,Vec X) 238 { 239 PetscInt i,j,Mx,My,xs,ys,xm,ym,its,l; 240 PetscReal hx,hy,hxdhy,hydhx; 241 PetscScalar **x,F,J,u,uxx,uyy; 242 DM da; 243 Vec localX; 244 245 PetscFunctionBeginUser; 246 PetscCall(SNESGetTolerances(snes,NULL,NULL,NULL,&its,NULL)); 247 PetscCall(SNESShellGetContext(snes,&da)); 248 249 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)); 250 251 hx = 1.0/(PetscReal)(Mx-1); 252 hy = 1.0/(PetscReal)(My-1); 253 hxdhy = hx/hy; 254 hydhx = hy/hx; 255 256 PetscCall(DMGetLocalVector(da,&localX)); 257 258 for (l=0; l<its; l++) { 259 260 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 261 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 262 /* 263 Get a pointer to vector data. 264 - For default PETSc vectors, VecGetArray() returns a pointer to 265 the data array. Otherwise, the routine is implementation dependent. 266 - You MUST call VecRestoreArray() when you no longer need access to 267 the array. 268 */ 269 PetscCall(DMDAVecGetArray(da,localX,&x)); 270 271 /* 272 Get local grid boundaries (for 2-dimensional DMDA): 273 xs, ys - starting grid indices (no ghost points) 274 xm, ym - widths of local grid (no ghost points) 275 276 */ 277 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 278 279 for (j=ys; j<ys+ym; j++) { 280 for (i=xs; i<xs+xm; i++) { 281 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 282 /* boundary conditions are all zero Dirichlet */ 283 x[j][i] = 0.0; 284 } else { 285 u = x[j][i]; 286 uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx; 287 uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy; 288 F = uxx + uyy; 289 J = 2.0*(hydhx + hxdhy); 290 u = u - F/J; 291 292 x[j][i] = u; 293 } 294 } 295 } 296 297 /* 298 Restore vector 299 */ 300 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 301 PetscCall(DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X)); 302 PetscCall(DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X)); 303 } 304 PetscCall(DMRestoreLocalVector(da,&localX)); 305 PetscFunctionReturn(0); 306 } 307 308 /*TEST 309 310 test: 311 args: -snes_monitor_short -snes_type nrichardson 312 requires: !single 313 314 test: 315 suffix: 2 316 args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale 317 requires: !single 318 319 test: 320 suffix: 3 321 args: -snes_monitor_short -snes_type ngmres 322 323 test: 324 suffix: 4 325 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none 326 327 test: 328 suffix: 5 329 args: -snes_monitor_short -snes_type ncg 330 331 test: 332 suffix: 6 333 args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none 334 335 test: 336 suffix: 7 337 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 338 requires: !single 339 340 test: 341 suffix: 8 342 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 343 344 test: 345 suffix: 9 346 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc 347 348 test: 349 suffix: 10 350 args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false 351 352 TEST*/ 353