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