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