1 2 static char help[] ="Model Equations for Advection \n"; 3 4 /* 5 Modified from ex3.c 6 Page 9, Section 1.2 Model Equations for Advection-Diffusion 7 8 u_t + a u_x = 0, 0<= x <= 1.0 9 10 The initial conditions used here different from the book. 11 12 Example: 13 ./ex6 -ts_monitor -ts_view_solution -ts_max_steps 100 -ts_monitor_solution draw -draw_pause .1 14 ./ex6 -ts_monitor -ts_max_steps 100 -ts_monitor_lg_error -draw_pause .1 15 */ 16 17 #include <petscts.h> 18 #include <petscdm.h> 19 #include <petscdmda.h> 20 21 /* 22 User-defined application context - contains data needed by the 23 application-provided call-back routines. 24 */ 25 typedef struct { 26 PetscReal a; /* advection strength */ 27 } AppCtx; 28 29 /* User-defined routines */ 30 extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*); 31 extern PetscErrorCode Solution(TS,PetscReal,Vec,AppCtx*); 32 extern PetscErrorCode IFunction_LaxFriedrichs(TS,PetscReal,Vec,Vec,Vec,void*); 33 extern PetscErrorCode IFunction_LaxWendroff(TS,PetscReal,Vec,Vec,Vec,void*); 34 35 int main(int argc,char **argv) 36 { 37 AppCtx appctx; /* user-defined application context */ 38 TS ts; /* timestepping context */ 39 Vec U; /* approximate solution vector */ 40 PetscErrorCode ierr; 41 PetscReal dt; 42 DM da; 43 PetscInt M; 44 PetscMPIInt rank; 45 PetscBool useLaxWendroff = PETSC_TRUE; 46 47 /* Initialize program and set problem parameters */ 48 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 49 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 50 51 appctx.a = -1.0; 52 ierr = PetscOptionsGetReal(NULL,NULL,"-a",&appctx.a,NULL);CHKERRQ(ierr); 53 54 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);CHKERRQ(ierr); 55 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 56 ierr = DMSetUp(da);CHKERRQ(ierr); 57 58 /* Create vector data structures for approximate and exact solutions */ 59 ierr = DMCreateGlobalVector(da,&U);CHKERRQ(ierr); 60 61 /* Create timestepping solver context */ 62 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 63 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 64 65 /* Function evaluation */ 66 ierr = PetscOptionsGetBool(NULL,NULL,"-useLaxWendroff",&useLaxWendroff,NULL);CHKERRQ(ierr); 67 if (useLaxWendroff) { 68 if (!rank) { 69 ierr = PetscPrintf(PETSC_COMM_SELF,"... Use Lax-Wendroff finite volume\n");CHKERRQ(ierr); 70 } 71 ierr = TSSetIFunction(ts,NULL,IFunction_LaxWendroff,&appctx);CHKERRQ(ierr); 72 } else { 73 if (!rank) { 74 ierr = PetscPrintf(PETSC_COMM_SELF,"... Use Lax-LaxFriedrichs finite difference\n");CHKERRQ(ierr); 75 } 76 ierr = TSSetIFunction(ts,NULL,IFunction_LaxFriedrichs,&appctx);CHKERRQ(ierr); 77 } 78 79 /* Customize timestepping solver */ 80 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 81 dt = 1.0/(PetscAbsReal(appctx.a)*M); 82 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 83 ierr = TSSetMaxSteps(ts,100);CHKERRQ(ierr); 84 ierr = TSSetMaxTime(ts,100.0);CHKERRQ(ierr); 85 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 86 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 87 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 88 89 /* Evaluate initial conditions */ 90 ierr = InitialConditions(ts,U,&appctx);CHKERRQ(ierr); 91 92 /* For testing accuracy of TS with already known solution, e.g., '-ts_monitor_lg_error' */ 93 ierr = TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);CHKERRQ(ierr); 94 95 /* Run the timestepping solver */ 96 ierr = TSSolve(ts,U);CHKERRQ(ierr); 97 98 /* Free work space */ 99 ierr = TSDestroy(&ts);CHKERRQ(ierr); 100 ierr = VecDestroy(&U);CHKERRQ(ierr); 101 ierr = DMDestroy(&da);CHKERRQ(ierr); 102 103 ierr = PetscFinalize(); 104 return ierr; 105 } 106 /* --------------------------------------------------------------------- */ 107 /* 108 InitialConditions - Computes the solution at the initial time. 109 110 Input Parameter: 111 u - uninitialized solution vector (global) 112 appctx - user-defined application context 113 114 Output Parameter: 115 u - vector with solution at initial time (global) 116 */ 117 PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx) 118 { 119 PetscScalar *u; 120 PetscErrorCode ierr; 121 PetscInt i,mstart,mend,um,M; 122 DM da; 123 PetscReal h; 124 125 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 126 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr); 127 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 128 h = 1.0/M; 129 mend = mstart + um; 130 /* 131 Get a pointer to vector data. 132 - For default PETSc vectors, VecGetArray() returns a pointer to 133 the data array. Otherwise, the routine is implementation dependent. 134 - You MUST call VecRestoreArray() when you no longer need access to 135 the array. 136 - Note that the Fortran interface to VecGetArray() differs from the 137 C version. See the users manual for details. 138 */ 139 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 140 141 /* 142 We initialize the solution array by simply writing the solution 143 directly into the array locations. Alternatively, we could use 144 VecSetValues() or VecSetValuesLocal(). 145 */ 146 for (i=mstart; i<mend; i++) u[i] = PetscSinReal(PETSC_PI*i*6.*h) + 3.*PetscSinReal(PETSC_PI*i*2.*h); 147 148 /* Restore vector */ 149 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 150 return 0; 151 } 152 /* --------------------------------------------------------------------- */ 153 /* 154 Solution - Computes the exact solution at a given time 155 156 Input Parameters: 157 t - current time 158 solution - vector in which exact solution will be computed 159 appctx - user-defined application context 160 161 Output Parameter: 162 solution - vector with the newly computed exact solution 163 u(x,t) = sin(6*PI*(x - a*t)) + 3 * sin(2*PI*(x - a*t)) 164 */ 165 PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx) 166 { 167 PetscScalar *u; 168 PetscReal a=appctx->a,h,PI6,PI2; 169 PetscErrorCode ierr; 170 PetscInt i,mstart,mend,um,M; 171 DM da; 172 173 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 174 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr); 175 ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 176 h = 1.0/M; 177 mend = mstart + um; 178 179 /* Get a pointer to vector data. */ 180 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 181 182 /* u[i] = sin(6*PI*(x[i] - a*t)) + 3 * sin(2*PI*(x[i] - a*t)) */ 183 PI6 = PETSC_PI*6.; 184 PI2 = PETSC_PI*2.; 185 for (i=mstart; i<mend; i++) { 186 u[i] = PetscSinReal(PI6*(i*h - a*t)) + 3.*PetscSinReal(PI2*(i*h - a*t)); 187 } 188 189 /* Restore vector */ 190 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 191 return 0; 192 } 193 194 /* --------------------------------------------------------------------- */ 195 /* 196 Use Lax-Friedrichs method to evaluate F(u,t) = du/dt + a * du/dx 197 198 See https://en.wikipedia.org/wiki/Lax%E2%80%93Friedrichs_method 199 */ 200 PetscErrorCode IFunction_LaxFriedrichs(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void* ctx) 201 { 202 PetscErrorCode ierr; 203 AppCtx *appctx=(AppCtx*)ctx; 204 PetscInt mstart,mend,M,i,um; 205 DM da; 206 Vec Uold,localUold; 207 PetscScalar *uarray,*f,*uoldarray,h,uave,c; 208 PetscReal dt; 209 210 PetscFunctionBegin; 211 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 212 ierr = TSGetSolution(ts,&Uold);CHKERRQ(ierr); 213 214 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 215 ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 216 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr); 217 h = 1.0/M; 218 mend = mstart + um; 219 /* printf(" mstart %d, um %d\n",mstart,um); */ 220 221 ierr = DMGetLocalVector(da,&localUold);CHKERRQ(ierr); 222 ierr = DMGlobalToLocalBegin(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr); 223 ierr = DMGlobalToLocalEnd(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr); 224 225 /* Get pointers to vector data */ 226 ierr = DMDAVecGetArrayRead(da,U,&uarray);CHKERRQ(ierr); 227 ierr = DMDAVecGetArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr); 228 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 229 230 /* advection */ 231 c = appctx->a*dt/h; /* Courant-Friedrichs-Lewy number (CFL number) */ 232 233 for (i=mstart; i<mend; i++) { 234 uave = 0.5*(uoldarray[i-1] + uoldarray[i+1]); 235 f[i] = uarray[i] - uave + c*0.5*(uoldarray[i+1] - uoldarray[i-1]); 236 } 237 238 /* Restore vectors */ 239 ierr = DMDAVecRestoreArrayRead(da,U,&uarray);CHKERRQ(ierr); 240 ierr = DMDAVecRestoreArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr); 241 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 242 ierr = DMRestoreLocalVector(da,&localUold);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 /* 247 Use Lax-Wendroff method to evaluate F(u,t) = du/dt + a * du/dx 248 */ 249 PetscErrorCode IFunction_LaxWendroff(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void* ctx) 250 { 251 PetscErrorCode ierr; 252 AppCtx *appctx=(AppCtx*)ctx; 253 PetscInt mstart,mend,M,i,um; 254 DM da; 255 Vec Uold,localUold; 256 PetscScalar *uarray,*f,*uoldarray,h,RFlux,LFlux,lambda; 257 PetscReal dt,a; 258 259 PetscFunctionBegin; 260 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 261 ierr = TSGetSolution(ts,&Uold);CHKERRQ(ierr); 262 263 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 264 ierr = DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 265 ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr); 266 h = 1.0/M; 267 mend = mstart + um; 268 /* printf(" mstart %d, um %d\n",mstart,um); */ 269 270 ierr = DMGetLocalVector(da,&localUold);CHKERRQ(ierr); 271 ierr = DMGlobalToLocalBegin(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr); 272 ierr = DMGlobalToLocalEnd(da,Uold,INSERT_VALUES,localUold);CHKERRQ(ierr); 273 274 /* Get pointers to vector data */ 275 ierr = DMDAVecGetArrayRead(da,U,&uarray);CHKERRQ(ierr); 276 ierr = DMDAVecGetArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr); 277 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 278 279 /* advection -- finite volume (appctx->a < 0 -- can be relaxed?) */ 280 lambda = dt/h; 281 a = appctx->a; 282 283 for (i=mstart; i<mend; i++) { 284 RFlux = 0.5 * a * (uoldarray[i+1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i+1] - uoldarray[i]); 285 LFlux = 0.5 * a * (uoldarray[i-1] + uoldarray[i]) - a*a*0.5*lambda * (uoldarray[i] - uoldarray[i-1]); 286 f[i] = uarray[i] - uoldarray[i] + lambda * (RFlux - LFlux); 287 } 288 289 /* Restore vectors */ 290 ierr = DMDAVecRestoreArrayRead(da,U,&uarray);CHKERRQ(ierr); 291 ierr = DMDAVecRestoreArrayRead(da,localUold,&uoldarray);CHKERRQ(ierr); 292 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 293 ierr = DMRestoreLocalVector(da,&localUold);CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 /*TEST 298 299 test: 300 args: -ts_max_steps 10 -ts_monitor 301 302 test: 303 suffix: 2 304 nsize: 3 305 args: -ts_max_steps 10 -ts_monitor 306 output_file: output/ex6_1.out 307 308 test: 309 suffix: 3 310 args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false 311 312 test: 313 suffix: 4 314 nsize: 3 315 args: -ts_max_steps 10 -ts_monitor -useLaxWendroff false 316 output_file: output/ex6_3.out 317 318 TEST*/ 319