1 static char help[] = "Surface processes in geophysics.\n\n"; 2 3 /*T 4 Concepts: SNES^parallel Surface process example 5 Concepts: DMDA^using distributed arrays; 6 Concepts: IS coloirng types; 7 Processors: n 8 T*/ 9 10 #include <petscsnes.h> 11 #include <petscdm.h> 12 #include <petscdmda.h> 13 14 /* 15 User-defined application context - contains data needed by the 16 application-provided call-back routines, FormJacobianLocal() and 17 FormFunctionLocal(). 18 */ 19 typedef struct { 20 PetscReal D; /* The diffusion coefficient */ 21 PetscReal K; /* The advection coefficient */ 22 PetscInt m; /* Exponent for A */ 23 } AppCtx; 24 25 /* 26 User-defined routines 27 */ 28 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*); 29 extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,AppCtx*); 30 31 int main(int argc,char **argv) 32 { 33 SNES snes; /* nonlinear solver */ 34 AppCtx user; /* user-defined work context */ 35 PetscInt its; /* iterations for convergence */ 36 PetscErrorCode ierr; 37 DM da; 38 39 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 40 Initialize program 41 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 42 43 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 44 45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 46 Initialize problem parameters 47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 48 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Surface Process Problem Options", "SNES");CHKERRQ(ierr); 49 user.D = 1.0; 50 CHKERRQ(PetscOptionsReal("-D", "The diffusion coefficient D", __FILE__, user.D, &user.D, NULL)); 51 user.K = 1.0; 52 CHKERRQ(PetscOptionsReal("-K", "The advection coefficient K", __FILE__, user.K, &user.K, NULL)); 53 user.m = 1; 54 CHKERRQ(PetscOptionsInt("-m", "The exponent for A", __FILE__, user.m, &user.m, NULL)); 55 ierr = PetscOptionsEnd();CHKERRQ(ierr); 56 57 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 58 Create distributed array (DMDA) to manage parallel grid and vectors 59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 60 CHKERRQ(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da)); 61 CHKERRQ(DMSetFromOptions(da)); 62 CHKERRQ(DMSetUp(da)); 63 CHKERRQ(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 64 CHKERRQ(DMSetApplicationContext(da,&user)); 65 CHKERRQ(SNESCreate(PETSC_COMM_WORLD, &snes)); 66 CHKERRQ(SNESSetDM(snes, da)); 67 68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 Set local function evaluation routine 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 CHKERRQ(DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user)); 72 73 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 Customize solver; set runtime options 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 CHKERRQ(SNESSetFromOptions(snes)); 77 78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 Solve nonlinear system 80 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 81 CHKERRQ(SNESSolve(snes,0,0)); 82 CHKERRQ(SNESGetIterationNumber(snes,&its)); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its)); 87 88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 89 Free work space. All PETSc objects should be destroyed when they 90 are no longer needed. 91 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 92 93 CHKERRQ(SNESDestroy(&snes)); 94 CHKERRQ(DMDestroy(&da)); 95 96 CHKERRQ(PetscFinalize()); 97 return 0; 98 } 99 100 PetscScalar funcU(DMDACoor2d *coords) 101 { 102 return coords->x + coords->y; 103 } 104 105 PetscScalar funcA(PetscScalar z, AppCtx *user) 106 { 107 PetscScalar v = 1.0; 108 PetscInt i; 109 110 for (i = 0; i < user->m; ++i) v *= z; 111 return v; 112 } 113 114 PetscScalar funcADer(PetscScalar z, AppCtx *user) 115 { 116 PetscScalar v = 1.0; 117 PetscInt i; 118 119 for (i = 0; i < user->m-1; ++i) v *= z; 120 return (PetscScalar)user->m*v; 121 } 122 123 /* 124 FormFunctionLocal - Evaluates nonlinear function, F(x). 125 */ 126 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user) 127 { 128 DM coordDA; 129 Vec coordinates; 130 DMDACoor2d **coords; 131 PetscScalar u, ux, uy, uxx, uyy; 132 PetscReal D, K, hx, hy, hxdhy, hydhx; 133 PetscInt i,j; 134 135 PetscFunctionBeginUser; 136 D = user->D; 137 K = user->K; 138 hx = 1.0/(PetscReal)(info->mx-1); 139 hy = 1.0/(PetscReal)(info->my-1); 140 hxdhy = hx/hy; 141 hydhx = hy/hx; 142 /* 143 Compute function over the locally owned part of the grid 144 */ 145 CHKERRQ(DMGetCoordinateDM(info->da, &coordDA)); 146 CHKERRQ(DMGetCoordinates(info->da, &coordinates)); 147 CHKERRQ(DMDAVecGetArray(coordDA, coordinates, &coords)); 148 for (j=info->ys; j<info->ys+info->ym; j++) { 149 for (i=info->xs; i<info->xs+info->xm; i++) { 150 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) f[j][i] = x[j][i]; 151 else { 152 u = x[j][i]; 153 ux = (x[j][i+1] - x[j][i])/hx; 154 uy = (x[j+1][i] - x[j][i])/hy; 155 uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx; 156 uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy; 157 f[j][i] = D*(uxx + uyy) - (K*funcA(x[j][i], user)*PetscSqrtScalar(ux*ux + uy*uy) + funcU(&coords[j][i]))*hx*hy; 158 PetscCheckFalse(PetscIsInfOrNanScalar(f[j][i]),PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", (double)PetscRealPart(f[j][i])); 159 } 160 } 161 } 162 CHKERRQ(DMDAVecRestoreArray(coordDA, coordinates, &coords)); 163 CHKERRQ(PetscLogFlops(11.0*info->ym*info->xm)); 164 PetscFunctionReturn(0); 165 } 166 167 /* 168 FormJacobianLocal - Evaluates Jacobian matrix. 169 */ 170 PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user) 171 { 172 MatStencil col[5], row; 173 PetscScalar D, K, A, v[5], hx, hy, hxdhy, hydhx, ux, uy; 174 PetscReal normGradZ; 175 PetscInt i, j,k; 176 177 PetscFunctionBeginUser; 178 D = user->D; 179 K = user->K; 180 hx = 1.0/(PetscReal)(info->mx-1); 181 hy = 1.0/(PetscReal)(info->my-1); 182 hxdhy = hx/hy; 183 hydhx = hy/hx; 184 185 /* 186 Compute entries for the locally owned part of the Jacobian. 187 - Currently, all PETSc parallel matrix formats are partitioned by 188 contiguous chunks of rows across the processors. 189 - Each processor needs to insert only elements that it owns 190 locally (but any non-local elements will be sent to the 191 appropriate processor during matrix assembly). 192 - Here, we set all entries for a particular row at once. 193 - We can set matrix entries either using either 194 MatSetValuesLocal() or MatSetValues(), as discussed above. 195 */ 196 for (j=info->ys; j<info->ys+info->ym; j++) { 197 for (i=info->xs; i<info->xs+info->xm; i++) { 198 row.j = j; row.i = i; 199 if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) { 200 /* boundary points */ 201 v[0] = 1.0; 202 CHKERRQ(MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES)); 203 } else { 204 /* interior grid points */ 205 ux = (x[j][i+1] - x[j][i])/hx; 206 uy = (x[j+1][i] - x[j][i])/hy; 207 normGradZ = PetscRealPart(PetscSqrtScalar(ux*ux + uy*uy)); 208 if (normGradZ < 1.0e-8) normGradZ = 1.0e-8; 209 A = funcA(x[j][i], user); 210 211 v[0] = -D*hxdhy; col[0].j = j - 1; col[0].i = i; 212 v[1] = -D*hydhx; col[1].j = j; col[1].i = i-1; 213 v[2] = D*2.0*(hydhx + hxdhy) + K*(funcADer(x[j][i], user)*normGradZ - A/normGradZ)*hx*hy; col[2].j = row.j; col[2].i = row.i; 214 v[3] = -D*hydhx + K*A*hx*hy/(2.0*normGradZ); col[3].j = j; col[3].i = i+1; 215 v[4] = -D*hxdhy + K*A*hx*hy/(2.0*normGradZ); col[4].j = j + 1; col[4].i = i; 216 for (k = 0; k < 5; ++k) { 217 PetscCheckFalse(PetscIsInfOrNanScalar(v[k]),PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", (double)PetscRealPart(v[k])); 218 } 219 CHKERRQ(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 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 229 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 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 CHKERRQ(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 235 PetscFunctionReturn(0); 236 } 237 238 /*TEST 239 240 test: 241 args: -snes_view -snes_monitor_short -da_refine 1 -pc_type mg -ksp_type fgmres -pc_mg_type full -mg_levels_ksp_chebyshev_esteig 0.5,1.1 242 243 test: 244 suffix: ew_1 245 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 1 246 requires: !single 247 248 test: 249 suffix: ew_2 250 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 2 251 252 test: 253 suffix: ew_3 254 args: -snes_monitor_short -ksp_converged_reason -da_grid_x 20 -da_grid_y 20 -snes_ksp_ew -snes_ksp_ew_version 3 255 requires: !single 256 257 test: 258 suffix: fm_rise_2 259 args: -K 3 -m 1 -D 0.2 -snes_monitor_short -snes_type ngmres -snes_npc_side right -npc_snes_type newtonls -npc_snes_linesearch_type basic -snes_ngmres_restart_it 1 -snes_ngmres_restart_fm_rise 260 requires: !single 261 262 test: 263 suffix: fm_rise_4 264 args: -K 3 -m 1 -D 0.2 -snes_monitor_short -snes_type ngmres -snes_npc_side right -npc_snes_type newtonls -npc_snes_linesearch_type basic -snes_ngmres_restart_it 2 -snes_ngmres_restart_fm_rise -snes_rtol 1.e-2 -snes_max_it 5 265 266 TEST*/ 267