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