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