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 ierr = PetscOptionsGetBool(NULL,NULL,"-usemonitor",&usemonitor,NULL);CHKERRQ(ierr); 49 50 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 51 Create distributed array (DMDA) to manage parallel grid and vectors 52 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 53 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr); 54 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 55 ierr = DMSetUp(da);CHKERRQ(ierr); 56 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr); 57 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr); 58 59 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60 Extract global vectors from DMDA; then duplicate for remaining 61 vectors that are the same types 62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 63 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 64 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 65 66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67 Create timestepping solver context 68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 69 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 70 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 71 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 72 ierr = TSSetRHSFunction(ts,NULL,FormFunction,da);CHKERRQ(ierr); 73 74 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 75 if (usemonitor) { 76 ierr = TSMonitorSet(ts,MyTSMonitor,0,0);CHKERRQ(ierr); 77 } 78 79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80 Customize nonlinear solver 81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 82 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 83 ierr = TSGetSNES(ts,&ts_snes);CHKERRQ(ierr); 84 if (usemonitor) { 85 ierr = PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);CHKERRQ(ierr); 86 ierr = SNESMonitorSet(ts_snes,(PetscErrorCode (*)(SNES,PetscInt,PetscReal,void *))MySNESMonitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);CHKERRQ(ierr); 87 } 88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 89 Set initial conditions 90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 91 ierr = FormInitialSolution(da,x);CHKERRQ(ierr); 92 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr); 93 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 94 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Set runtime options 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Solve nonlinear system 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 ierr = TSSolve(ts,x);CHKERRQ(ierr); 105 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 106 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 107 ierr = VecViewFromOptions(x,NULL,"-final_sol");CHKERRQ(ierr); 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Free work space. All PETSc objects should be destroyed when they 111 are no longer needed. 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 ierr = VecDestroy(&x);CHKERRQ(ierr); 114 ierr = VecDestroy(&r);CHKERRQ(ierr); 115 ierr = TSDestroy(&ts);CHKERRQ(ierr); 116 ierr = DMDestroy(&da);CHKERRQ(ierr); 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 PetscErrorCode ierr; 137 PetscInt i,j,Mx,My,xs,ys,xm,ym; 138 PetscReal hx,hy,/*hxdhy,hydhx,*/ sx,sy; 139 PetscScalar u,uxx,uyy,v,***x,***f; 140 Vec localX; 141 142 PetscFunctionBeginUser; 143 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 144 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); 145 146 hx = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx); 147 hy = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy); 148 /*hxdhy = hx/hy;*/ 149 /*hydhx = hy/hx;*/ 150 151 /* 152 Scatter ghost points to local vector,using the 2-step process 153 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 154 By placing code between these two statements, computations can be 155 done while messages are in transition. 156 */ 157 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 158 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 159 160 /* 161 Get pointers to vector data 162 */ 163 ierr = DMDAVecGetArrayDOF(da,localX,&x);CHKERRQ(ierr); 164 ierr = DMDAVecGetArrayDOF(da,F,&f);CHKERRQ(ierr); 165 166 /* 167 Get local grid boundaries 168 */ 169 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 170 171 /* 172 Compute function over the locally owned part of the grid 173 */ 174 for (j=ys; j<ys+ym; j++) { 175 for (i=xs; i<xs+xm; i++) { 176 if (i == 0 || j == 0 || i == Mx-1 || j == My-1) { 177 f[j][i][0] = x[j][i][0]; 178 f[j][i][1] = x[j][i][1]; 179 continue; 180 } 181 u = x[j][i][0]; 182 v = x[j][i][1]; 183 uxx = (-2.0*u + x[j][i-1][0] + x[j][i+1][0])*sx; 184 uyy = (-2.0*u + x[j-1][i][0] + x[j+1][i][0])*sy; 185 f[j][i][0] = v; 186 f[j][i][1] = uxx + uyy; 187 } 188 } 189 190 /* 191 Restore vectors 192 */ 193 ierr = DMDAVecRestoreArrayDOF(da,localX,&x);CHKERRQ(ierr); 194 ierr = DMDAVecRestoreArrayDOF(da,F,&f);CHKERRQ(ierr); 195 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 196 ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 /* ------------------------------------------------------------------- */ 201 PetscErrorCode FormInitialSolution(DM da,Vec U) 202 { 203 PetscErrorCode ierr; 204 PetscInt i,j,xs,ys,xm,ym,Mx,My; 205 PetscScalar ***u; 206 PetscReal hx,hy,x,y,r; 207 208 PetscFunctionBeginUser; 209 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); 210 211 hx = 1.0/(PetscReal)(Mx-1); 212 hy = 1.0/(PetscReal)(My-1); 213 214 /* 215 Get pointers to vector data 216 */ 217 ierr = DMDAVecGetArrayDOF(da,U,&u);CHKERRQ(ierr); 218 219 /* 220 Get local grid boundaries 221 */ 222 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 223 224 /* 225 Compute function over the locally owned part of the grid 226 */ 227 for (j=ys; j<ys+ym; j++) { 228 y = j*hy; 229 for (i=xs; i<xs+xm; i++) { 230 x = i*hx; 231 r = PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)); 232 if (r < .125) { 233 u[j][i][0] = PetscExpReal(-30.0*r*r*r); 234 u[j][i][1] = 0.0; 235 } else { 236 u[j][i][0] = 0.0; 237 u[j][i][1] = 0.0; 238 } 239 } 240 } 241 242 /* 243 Restore vectors 244 */ 245 ierr = DMDAVecRestoreArrayDOF(da,U,&u);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 PetscErrorCode MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,void *ctx) 250 { 251 PetscErrorCode ierr; 252 PetscReal norm; 253 MPI_Comm comm; 254 255 PetscFunctionBeginUser; 256 ierr = VecNorm(v,NORM_2,&norm);CHKERRQ(ierr); 257 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 258 if (step > -1) { /* -1 is used to indicate an interpolated value */ 259 ierr = PetscPrintf(comm,"timestep %D time %g norm %g\n",step,(double)ptime,(double)norm);CHKERRQ(ierr); 260 } 261 PetscFunctionReturn(0); 262 } 263 264 /* 265 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES. 266 Input Parameters: 267 snes - the SNES context 268 its - iteration number 269 fnorm - 2-norm function value (may be estimated) 270 ctx - optional user-defined context for private data for the 271 monitor routine, as set by SNESMonitorSet() 272 */ 273 PetscErrorCode MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *vf) 274 { 275 PetscErrorCode ierr; 276 277 PetscFunctionBeginUser; 278 ierr = SNESMonitorDefaultShort(snes,its,fnorm,vf);CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 /*TEST 282 283 test: 284 args: -da_grid_x 20 -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short 285 requires: !single 286 287 test: 288 suffix: 2 289 args: -da_grid_x 20 -ts_max_time 0.11 -ts_dt 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short 290 requires: !single 291 292 test: 293 suffix: glvis_da_2d_vect 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 295 requires: !single 296 297 test: 298 suffix: glvis_da_2d_vect_ll 299 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 300 requires: !single 301 302 TEST*/ 303