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