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 DM da; 31 PetscViewerAndFormat *vf; 32 33 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 34 Initialize program 35 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 36 PetscFunctionBeginUser; 37 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 38 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 39 Create distributed array (DMDA) to manage parallel grid and vectors 40 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 41 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da)); 42 PetscCall(DMSetFromOptions(da)); 43 PetscCall(DMSetUp(da)); 44 45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 46 Extract global vectors from DMDA; then duplicate for remaining 47 vectors that are the same types 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 PetscCall(DMCreateGlobalVector(da,&x)); 50 PetscCall(VecDuplicate(x,&r)); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create timestepping solver context 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 56 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 57 PetscCall(TSSetRHSFunction(ts,NULL,FormFunction,da)); 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 PetscCall(TSSetMaxTime(ts,1.0)); 73 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 74 PetscCall(TSMonitorSet(ts,MyTSMonitor,PETSC_VIEWER_STDOUT_WORLD,NULL)); 75 PetscCall(TSSetDM(ts,da)); 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Customize nonlinear solver 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 79 PetscCall(TSSetType(ts,TSBEULER)); 80 PetscCall(TSGetSNES(ts,&snes)); 81 PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf)); 82 PetscCall(SNESMonitorSet(snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void*))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy)); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Set initial conditions 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 PetscCall(FormInitialSolution(da,x)); 88 PetscCall(TSSetTimeStep(ts,.0001)); 89 PetscCall(TSSetSolution(ts,x)); 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Set runtime options 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 PetscCall(TSSetFromOptions(ts)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Solve nonlinear system 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 PetscCall(TSSolve(ts,x)); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Free work space. All PETSc objects should be destroyed when they 103 are no longer needed. 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 PetscCall(VecDestroy(&x)); 106 PetscCall(VecDestroy(&r)); 107 PetscCall(TSDestroy(&ts)); 108 PetscCall(DMDestroy(&da)); 109 110 PetscCall(PetscFinalize()); 111 return 0; 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 PetscInt i,j,Mx,My,xs,ys,xm,ym; 129 PetscReal two = 2.0,hx,hy,sx,sy; 130 PetscScalar u,uxx,uyy,**x,**f; 131 Vec localX; 132 133 PetscFunctionBeginUser; 134 PetscCall(TSGetDM(ts,&da)); 135 PetscCall(DMGetLocalVector(da,&localX)); 136 PetscCall(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)); 137 138 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 139 hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy); 140 141 /* 142 Scatter ghost points to local vector,using the 2-step process 143 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 144 By placing code between these two statements, computations can be 145 done while messages are in transition. 146 */ 147 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 148 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 149 150 /* 151 Get pointers to vector data 152 */ 153 PetscCall(DMDAVecGetArrayRead(da,localX,&x)); 154 PetscCall(DMDAVecGetArray(da,F,&f)); 155 156 /* 157 Get local grid boundaries 158 */ 159 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 160 161 /* 162 Compute function over the locally owned part of the grid 163 */ 164 for (j=ys; j<ys+ym; j++) { 165 for (i=xs; i<xs+xm; i++) { 166 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 167 f[j][i] = x[j][i]; 168 continue; 169 } 170 u = x[j][i]; 171 uxx = (two*u - x[j][i-1] - x[j][i+1])*sx; 172 uyy = (two*u - x[j-1][i] - x[j+1][i])*sy; 173 /* f[j][i] = -(uxx + uyy); */ 174 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 + 175 (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy); 176 } 177 } 178 179 /* 180 Restore vectors 181 */ 182 PetscCall(DMDAVecRestoreArrayRead(da,localX,&x)); 183 PetscCall(DMDAVecRestoreArray(da,F,&f)); 184 PetscCall(DMRestoreLocalVector(da,&localX)); 185 PetscCall(PetscLogFlops(11.0*ym*xm)); 186 PetscFunctionReturn(0); 187 } 188 189 /* ------------------------------------------------------------------- */ 190 PetscErrorCode FormInitialSolution(DM da,Vec U) 191 { 192 PetscInt i,j,xs,ys,xm,ym,Mx,My; 193 PetscScalar **u; 194 PetscReal hx,hy,x,y,r; 195 196 PetscFunctionBeginUser; 197 PetscCall(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)); 198 199 hx = 1.0/(PetscReal)(Mx-1); 200 hy = 1.0/(PetscReal)(My-1); 201 202 /* 203 Get pointers to vector data 204 */ 205 PetscCall(DMDAVecGetArray(da,U,&u)); 206 207 /* 208 Get local grid boundaries 209 */ 210 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 211 212 /* 213 Compute function over the locally owned part of the grid 214 */ 215 for (j=ys; j<ys+ym; j++) { 216 y = j*hy; 217 for (i=xs; i<xs+xm; i++) { 218 x = i*hx; 219 r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)); 220 if (r < .125) u[j][i] = PetscExpReal(-30.0*r*r*r); 221 else u[j][i] = 0.0; 222 } 223 } 224 225 /* 226 Restore vectors 227 */ 228 PetscCall(DMDAVecRestoreArray(da,U,&u)); 229 PetscFunctionReturn(0); 230 } 231 232 PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx) 233 { 234 PetscReal norm; 235 MPI_Comm comm; 236 237 PetscFunctionBeginUser; 238 if (step < 0) PetscFunctionReturn(0); /* step of -1 indicates an interpolated solution */ 239 PetscCall(VecNorm(v,NORM_2,&norm)); 240 PetscCall(PetscObjectGetComm((PetscObject)ts,&comm)); 241 PetscCall(PetscPrintf(comm,"timestep %" PetscInt_FMT " time %g norm %g\n",step,(double)ptime,(double)norm)); 242 PetscFunctionReturn(0); 243 } 244 245 /* 246 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES. 247 Input Parameters: 248 snes - the SNES context 249 its - iteration number 250 fnorm - 2-norm function value (may be estimated) 251 ctx - optional user-defined context for private data for the 252 monitor routine, as set by SNESMonitorSet() 253 */ 254 PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf) 255 { 256 PetscFunctionBeginUser; 257 PetscCall(SNESMonitorDefaultShort(snes,its,fnorm,vf)); 258 PetscFunctionReturn(0); 259 } 260 261 /*TEST 262 263 test: 264 args: -ts_max_steps 5 265 266 test: 267 suffix: 2 268 args: -ts_max_steps 5 -snes_mf_operator 269 270 test: 271 suffix: 3 272 args: -ts_max_steps 5 -snes_mf -pc_type none 273 274 TEST*/ 275