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