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