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