1 2 static char help[] = "Nonlinear, time-dependent PDE in 2d.\n"; 3 4 /* 5 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 6 Include "petscts.h" so that we can use SNES solvers. Note that this 7 file automatically includes: 8 petscsys.h - base PETSc routines petscvec.h - vectors 9 petscmat.h - matrices 10 petscis.h - index sets petscksp.h - Krylov subspace methods 11 petscviewer.h - viewers petscpc.h - preconditioners 12 petscksp.h - linear solvers 13 */ 14 #include <petscdm.h> 15 #include <petscdmda.h> 16 #include <petscts.h> 17 18 /* 19 User-defined routines 20 */ 21 extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec); 22 extern PetscErrorCode MyTSMonitor(TS,PetscInt,PetscReal,Vec,void*); 23 extern PetscErrorCode MySNESMonitor(SNES,PetscInt,PetscReal,PetscViewerAndFormat*); 24 25 int main(int argc,char **argv) 26 { 27 TS ts; /* time integrator */ 28 SNES snes; 29 Vec x,r; /* solution, residual vectors */ 30 PetscErrorCode ierr; 31 DM da; 32 PetscViewerAndFormat *vf; 33 34 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35 Initialize program 36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 37 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 38 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 39 Create distributed array (DMDA) to manage parallel grid and vectors 40 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 41 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr); 42 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 43 ierr = DMSetUp(da);CHKERRQ(ierr); 44 45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 46 Extract global vectors from DMDA; then duplicate for remaining 47 vectors that are the same types 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 50 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create timestepping solver context 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 56 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 57 ierr = TSSetRHSFunction(ts,NULL,FormFunction,da);CHKERRQ(ierr); 58 59 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60 Create matrix data structure; set Jacobian evaluation routine 61 62 Set Jacobian matrix data structure and default Jacobian evaluation 63 routine. User can override with: 64 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 65 (unless user explicitly sets preconditioner) 66 -snes_mf_operator : form preconditioning matrix as set by the user, 67 but use matrix-free approx for Jacobian-vector 68 products within Newton-Krylov method 69 70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 71 72 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 73 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 74 ierr = TSMonitorSet(ts,MyTSMonitor,PETSC_VIEWER_STDOUT_WORLD,NULL);CHKERRQ(ierr); 75 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Customize nonlinear solver 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 79 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 80 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 81 ierr = PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);CHKERRQ(ierr); 82 ierr = SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);CHKERRQ(ierr); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Set initial conditions 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 ierr = FormInitialSolution(da,x);CHKERRQ(ierr); 88 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr); 89 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Set runtime options 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Solve nonlinear system 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 ierr = TSSolve(ts,x);CHKERRQ(ierr); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Free work space. All PETSc objects should be destroyed when they 103 are no longer needed. 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 ierr = VecDestroy(&x);CHKERRQ(ierr); 106 ierr = VecDestroy(&r);CHKERRQ(ierr); 107 ierr = TSDestroy(&ts);CHKERRQ(ierr); 108 ierr = DMDestroy(&da);CHKERRQ(ierr); 109 110 ierr = PetscFinalize(); 111 return ierr; 112 } 113 /* ------------------------------------------------------------------- */ 114 /* 115 FormFunction - Evaluates nonlinear function, F(x). 116 117 Input Parameters: 118 . ts - the TS context 119 . X - input vector 120 . ptr - optional user-defined context, as set by SNESSetFunction() 121 122 Output Parameter: 123 . F - function vector 124 */ 125 PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr) 126 { 127 DM da; 128 PetscErrorCode ierr; 129 PetscInt i,j,Mx,My,xs,ys,xm,ym; 130 PetscReal two = 2.0,hx,hy,sx,sy; 131 PetscScalar u,uxx,uyy,**x,**f; 132 Vec localX; 133 134 PetscFunctionBeginUser; 135 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 136 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 137 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); 138 139 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 140 hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy); 141 142 /* 143 Scatter ghost points to local vector,using the 2-step process 144 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 145 By placing code between these two statements, computations can be 146 done while messages are in transition. 147 */ 148 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 149 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 150 151 /* 152 Get pointers to vector data 153 */ 154 ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); 155 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 156 157 /* 158 Get local grid boundaries 159 */ 160 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 161 162 /* 163 Compute function over the locally owned part of the grid 164 */ 165 for (j=ys; j<ys+ym; j++) { 166 for (i=xs; i<xs+xm; i++) { 167 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 168 f[j][i] = x[j][i]; 169 continue; 170 } 171 u = x[j][i]; 172 uxx = (two*u - x[j][i-1] - x[j][i+1])*sx; 173 uyy = (two*u - x[j-1][i] - x[j+1][i])*sy; 174 /* f[j][i] = -(uxx + uyy); */ 175 f[j][i] = -u*(uxx + uyy) - (4.0 - 1.0)*((x[j][i+1] - x[j][i-1])*(x[j][i+1] - x[j][i-1])*.25*sx + 176 (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy); 177 } 178 } 179 180 /* 181 Restore vectors 182 */ 183 ierr = DMDAVecRestoreArrayRead(da,localX,&x);CHKERRQ(ierr); 184 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 185 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 186 ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 /* ------------------------------------------------------------------- */ 191 PetscErrorCode FormInitialSolution(DM da,Vec U) 192 { 193 PetscErrorCode ierr; 194 PetscInt i,j,xs,ys,xm,ym,Mx,My; 195 PetscScalar **u; 196 PetscReal hx,hy,x,y,r; 197 198 PetscFunctionBeginUser; 199 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); 200 201 hx = 1.0/(PetscReal)(Mx-1); 202 hy = 1.0/(PetscReal)(My-1); 203 204 /* 205 Get pointers to vector data 206 */ 207 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 208 209 /* 210 Get local grid boundaries 211 */ 212 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 213 214 /* 215 Compute function over the locally owned part of the grid 216 */ 217 for (j=ys; j<ys+ym; j++) { 218 y = j*hy; 219 for (i=xs; i<xs+xm; i++) { 220 x = i*hx; 221 r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)); 222 if (r < .125) u[j][i] = PetscExpReal(-30.0*r*r*r); 223 else u[j][i] = 0.0; 224 } 225 } 226 227 /* 228 Restore vectors 229 */ 230 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx) 235 { 236 PetscErrorCode ierr; 237 PetscReal norm; 238 MPI_Comm comm; 239 240 PetscFunctionBeginUser; 241 if (step < 0) PetscFunctionReturn(0); /* step of -1 indicates an interpolated solution */ 242 ierr = VecNorm(v,NORM_2,&norm);CHKERRQ(ierr); 243 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 244 ierr = PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);CHKERRQ(ierr); 245 PetscFunctionReturn(0); 246 } 247 248 /* 249 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES. 250 Input Parameters: 251 snes - the SNES context 252 its - iteration number 253 fnorm - 2-norm function value (may be estimated) 254 ctx - optional user-defined context for private data for the 255 monitor routine, as set by SNESMonitorSet() 256 */ 257 PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf) 258 { 259 PetscErrorCode ierr; 260 261 PetscFunctionBeginUser; 262 ierr = SNESMonitorDefaultShort(snes,its,fnorm,vf);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 /*TEST 267 268 test: 269 args: -ts_max_steps 5 270 271 test: 272 suffix: 2 273 args: -ts_max_steps 5 -snes_mf_operator 274 275 test: 276 suffix: 3 277 args: -ts_max_steps 5 -snes_mf -pc_type none 278 279 TEST*/ 280