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