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