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