1 2 3 static char help[] ="Model Equations for Advection-Diffusion\n"; 4 5 /* 6 Page 9, Section 1.2 Model Equations for Advection-Diffusion 7 8 u_t = a u_x + d u_xx 9 10 The initial conditions used here different then in the book. 11 12 */ 13 14 /* 15 Helpful runtime linear solver options: 16 -pc_type mg -da_refine 2 -snes_monitor -ksp_monitor -ts_view (geometric multigrid with three levels) 17 18 */ 19 20 /* 21 Include "petscts.h" so that we can use TS solvers. Note that this file 22 automatically includes: 23 petscsys.h - base PETSc routines petscvec.h - vectors 24 petscmat.h - matrices 25 petscis.h - index sets petscksp.h - Krylov subspace methods 26 petscviewer.h - viewers petscpc.h - preconditioners 27 petscksp.h - linear solvers petscsnes.h - nonlinear solvers 28 */ 29 30 #include <petscts.h> 31 #include <petscdm.h> 32 #include <petscdmda.h> 33 34 /* 35 User-defined application context - contains data needed by the 36 application-provided call-back routines. 37 */ 38 typedef struct { 39 PetscScalar a,d; /* advection and diffusion strength */ 40 PetscBool upwind; 41 } AppCtx; 42 43 /* 44 User-defined routines 45 */ 46 extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*); 47 extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*); 48 extern PetscErrorCode Solution(TS,PetscReal,Vec,AppCtx*); 49 50 int main(int argc,char **argv) 51 { 52 AppCtx appctx; /* user-defined application context */ 53 TS ts; /* timestepping context */ 54 Vec U; /* approximate solution vector */ 55 PetscErrorCode ierr; 56 PetscReal dt; 57 DM da; 58 PetscInt M; 59 60 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 61 Initialize program and set problem parameters 62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 63 64 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 65 appctx.a = 1.0; 66 appctx.d = 0.0; 67 ierr = PetscOptionsGetScalar(NULL,NULL,"-a",&appctx.a,NULL);CHKERRQ(ierr); 68 ierr = PetscOptionsGetScalar(NULL,NULL,"-d",&appctx.d,NULL);CHKERRQ(ierr); 69 appctx.upwind = PETSC_TRUE; 70 ierr = PetscOptionsGetBool(NULL,NULL,"-upwind",&appctx.upwind,NULL);CHKERRQ(ierr); 71 72 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);CHKERRQ(ierr); 73 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 74 ierr = DMSetUp(da);CHKERRQ(ierr); 75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 Create vector data structures 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 79 /* 80 Create vector data structures for approximate and exact solutions 81 */ 82 ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Create timestepping solver context 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 87 88 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 89 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 90 91 /* 92 For linear problems with a time-dependent f(U,t) in the equation 93 u_t = f(u,t), the user provides the discretized right-hand-side 94 as a time-dependent matrix. 95 */ 96 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);CHKERRQ(ierr); 97 ierr = TSSetRHSJacobian(ts,NULL,NULL,RHSMatrixHeat,&appctx);CHKERRQ(ierr); 98 ierr = TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);CHKERRQ(ierr); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Customize timestepping solver: 102 - Set timestepping duration info 103 Then set runtime options, which can override these defaults. 104 For example, 105 -ts_max_steps <maxsteps> -ts_max_time <maxtime> 106 to override the defaults set by TSSetMaxSteps()/TSSetMaxTime(). 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 109 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 110 dt = .48/(M*M); 111 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 112 ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 113 ierr = TSSetMaxTime(ts,100.0);CHKERRQ(ierr); 114 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 115 ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); 116 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 117 118 /* 119 Evaluate initial conditions 120 */ 121 ierr = InitialConditions(ts,U,&appctx);CHKERRQ(ierr); 122 123 /* 124 Run the timestepping solver 125 */ 126 ierr = TSSolve(ts,U);CHKERRQ(ierr); 127 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Free work space. All PETSc objects should be destroyed when they 131 are no longer needed. 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 134 ierr = TSDestroy(&ts);CHKERRQ(ierr); 135 ierr = VecDestroy(&U);CHKERRQ(ierr); 136 ierr = DMDestroy(&da);CHKERRQ(ierr); 137 138 /* 139 Always call PetscFinalize() before exiting a program. This routine 140 - finalizes the PETSc libraries as well as MPI 141 - provides summary and diagnostic information if certain runtime 142 options are chosen (e.g., -log_view). 143 */ 144 ierr = PetscFinalize(); 145 return ierr; 146 } 147 /* --------------------------------------------------------------------- */ 148 /* 149 InitialConditions - Computes the solution at the initial time. 150 151 Input Parameter: 152 u - uninitialized solution vector (global) 153 appctx - user-defined application context 154 155 Output Parameter: 156 u - vector with solution at initial time (global) 157 */ 158 PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx) 159 { 160 PetscScalar *u,h; 161 PetscErrorCode ierr; 162 PetscInt i,mstart,mend,xm,M; 163 DM da; 164 165 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 166 ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr); 167 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 168 h = 1.0/M; 169 mend = mstart + xm; 170 /* 171 Get a pointer to vector data. 172 - For default PETSc vectors, VecGetArray() returns a pointer to 173 the data array. Otherwise, the routine is implementation dependent. 174 - You MUST call VecRestoreArray() when you no longer need access to 175 the array. 176 - Note that the Fortran interface to VecGetArray() differs from the 177 C version. See the users manual for details. 178 */ 179 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 180 181 /* 182 We initialize the solution array by simply writing the solution 183 directly into the array locations. Alternatively, we could use 184 VecSetValues() or VecSetValuesLocal(). 185 */ 186 for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h); 187 188 /* 189 Restore vector 190 */ 191 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 192 return 0; 193 } 194 /* --------------------------------------------------------------------- */ 195 /* 196 Solution - Computes the exact solution at a given time. 197 198 Input Parameters: 199 t - current time 200 solution - vector in which exact solution will be computed 201 appctx - user-defined application context 202 203 Output Parameter: 204 solution - vector with the newly computed exact solution 205 */ 206 PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx) 207 { 208 PetscScalar *u,ex1,ex2,sc1,sc2,h; 209 PetscErrorCode ierr; 210 PetscInt i,mstart,mend,xm,M; 211 DM da; 212 213 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 214 ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr); 215 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 216 h = 1.0/M; 217 mend = mstart + xm; 218 /* 219 Get a pointer to vector data. 220 */ 221 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 222 223 /* 224 Simply write the solution directly into the array locations. 225 Alternatively, we culd use VecSetValues() or VecSetValuesLocal(). 226 */ 227 ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t); 228 ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t); 229 sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h; 230 for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(sc1*(PetscReal)i + appctx->a*PETSC_PI*6.*t)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i + appctx->a*PETSC_PI*2.*t)*ex2; 231 232 /* 233 Restore vector 234 */ 235 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 236 return 0; 237 } 238 239 /* --------------------------------------------------------------------- */ 240 /* 241 RHSMatrixHeat - User-provided routine to compute the right-hand-side 242 matrix for the heat equation. 243 244 Input Parameters: 245 ts - the TS context 246 t - current time 247 global_in - global input vector 248 dummy - optional user-defined context, as set by TSetRHSJacobian() 249 250 Output Parameters: 251 AA - Jacobian matrix 252 BB - optionally different preconditioning matrix 253 str - flag indicating matrix structure 254 255 Notes: 256 Recall that MatSetValues() uses 0-based row and column numbers 257 in Fortran as well as in C. 258 */ 259 PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat AA,Mat BB,void *ctx) 260 { 261 Mat A = AA; /* Jacobian matrix */ 262 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 263 PetscInt mstart, mend; 264 PetscErrorCode ierr; 265 PetscInt i,idx[3],M,xm; 266 PetscScalar v[3],h; 267 DM da; 268 269 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 270 ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 271 ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr); 272 h = 1.0/M; 273 mend = mstart + xm; 274 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 275 Compute entries for the locally owned part of the matrix 276 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 277 /* 278 Set matrix rows corresponding to boundary data 279 */ 280 281 /* diffusion */ 282 v[0] = appctx->d/(h*h); 283 v[1] = -2.0*appctx->d/(h*h); 284 v[2] = appctx->d/(h*h); 285 if (!mstart) { 286 idx[0] = M-1; idx[1] = 0; idx[2] = 1; 287 ierr = MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 288 mstart++; 289 } 290 291 if (mend == M) { 292 mend--; 293 idx[0] = M-2; idx[1] = M-1; idx[2] = 0; 294 ierr = MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 295 } 296 297 /* 298 Set matrix rows corresponding to interior data. We construct the 299 matrix one row at a time. 300 */ 301 for (i=mstart; i<mend; i++) { 302 idx[0] = i-1; idx[1] = i; idx[2] = i+1; 303 ierr = MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 304 } 305 ierr = MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 306 ierr = MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr); 307 308 ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr); 309 mend = mstart + xm; 310 if (!appctx->upwind) { 311 /* advection -- centered differencing */ 312 v[0] = -.5*appctx->a/(h); 313 v[1] = .5*appctx->a/(h); 314 if (!mstart) { 315 idx[0] = M-1; idx[1] = 1; 316 ierr = MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 317 mstart++; 318 } 319 320 if (mend == M) { 321 mend--; 322 idx[0] = M-2; idx[1] = 0; 323 ierr = MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 324 } 325 326 for (i=mstart; i<mend; i++) { 327 idx[0] = i-1; idx[1] = i+1; 328 ierr = MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 329 } 330 } else { 331 /* advection -- upwinding */ 332 v[0] = -appctx->a/(h); 333 v[1] = appctx->a/(h); 334 if (!mstart) { 335 idx[0] = 0; idx[1] = 1; 336 ierr = MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 337 mstart++; 338 } 339 340 if (mend == M) { 341 mend--; 342 idx[0] = M-1; idx[1] = 0; 343 ierr = MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 344 } 345 346 for (i=mstart; i<mend; i++) { 347 idx[0] = i; idx[1] = i+1; 348 ierr = MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);CHKERRQ(ierr); 349 } 350 } 351 352 353 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 354 Complete the matrix assembly process and set some options 355 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 356 /* 357 Assemble matrix, using the 2-step process: 358 MatAssemblyBegin(), MatAssemblyEnd() 359 Computations can be done while messages are in transition 360 by placing code between these two statements. 361 */ 362 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 363 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 364 365 /* 366 Set and option to indicate that we will never add a new nonzero location 367 to the matrix. If we do, it will generate an error. 368 */ 369 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 370 return 0; 371 } 372 373 374 /*TEST 375 376 test: 377 args: -pc_type mg -da_refine 2 -ts_view -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3 378 requires: double 379 filter: grep -v "total number of" 380 381 test: 382 suffix: 2 383 args: -pc_type mg -da_refine 2 -ts_view -ts_monitor_draw_solution -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3 384 requires: x 385 output_file: output/ex3_1.out 386 requires: double 387 filter: grep -v "total number of" 388 389 TEST*/ 390