1 /* 2 Note: 3 -hratio is the ratio between mesh size of corse grids and fine grids 4 */ 5 6 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 7 " advect - Constant coefficient scalar advection\n" 8 " u_t + (a*u)_x = 0\n" 9 " shallow - 1D Shallow water equations (Saint Venant System)\n" 10 " h_t + (q)_x = 0 \n" 11 " q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0 \n" 12 " where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n" 13 " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n" 14 " hxs = hratio*hxf \n" 15 " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n" 16 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n" 17 " the states across shocks and rarefactions\n" 18 " simulation - use reference solution which is generated by smaller time step size to be true solution,\n" 19 " also the reference solution should be generated by user and stored in a binary file.\n" 20 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n" 21 " bc_type - Boundary condition for the problem, options are: periodic, outflow, inflow " 22 "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n" 23 "The problem size should be set with -da_grid_x M\n\n"; 24 25 /* 26 Example: 27 ./ex4 -da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 28 ./ex4 -da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 2.5 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 29 ./ex4 -da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 30 ./ex4 -da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_dt 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 31 ./ex4 -da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_dt 0.01 -ts_max_time 5.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 32 33 Contributed by: Aidan Hamilton <aidan@udel.edu> 34 */ 35 36 #include <petscts.h> 37 #include <petscdm.h> 38 #include <petscdmda.h> 39 #include <petscdraw.h> 40 #include "finitevolume1d.h" 41 #include <petsc/private/kernels/blockinvert.h> 42 43 static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); } 44 static inline PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; } 45 /* --------------------------------- Advection ----------------------------------- */ 46 typedef struct { 47 PetscReal a; /* advective velocity */ 48 } AdvectCtx; 49 50 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 51 { 52 AdvectCtx *ctx = (AdvectCtx*)vctx; 53 PetscReal speed; 54 55 PetscFunctionBeginUser; 56 speed = ctx->a; 57 flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0]; 58 *maxspeed = speed; 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 63 { 64 AdvectCtx *ctx = (AdvectCtx*)vctx; 65 66 PetscFunctionBeginUser; 67 X[0] = 1.; 68 Xi[0] = 1.; 69 speeds[0] = ctx->a; 70 PetscFunctionReturn(0); 71 } 72 73 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 74 { 75 AdvectCtx *ctx = (AdvectCtx*)vctx; 76 PetscReal a = ctx->a,x0; 77 78 PetscFunctionBeginUser; 79 switch (bctype) { 80 case FVBC_OUTFLOW: x0 = x-a*t; break; 81 case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break; 82 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 83 } 84 switch (initial) { 85 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 86 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 87 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 88 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break; 89 case 4: u[0] = PetscAbs(x0); break; 90 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break; 91 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break; 92 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break; 93 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 94 } 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 99 { 100 AdvectCtx *user; 101 102 PetscFunctionBeginUser; 103 PetscCall(PetscNew(&user)); 104 ctx->physics2.sample2 = PhysicsSample_Advect; 105 ctx->physics2.riemann2 = PhysicsRiemann_Advect; 106 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; 107 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; 108 ctx->physics2.user = user; 109 ctx->physics2.dof = 1; 110 111 PetscCall(PetscStrallocpy("u",&ctx->physics2.fieldname[0])); 112 user->a = 1; 113 PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection",""); 114 PetscCall(PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL)); 115 PetscOptionsEnd(); 116 PetscFunctionReturn(0); 117 } 118 /* --------------------------------- Shallow Water ----------------------------------- */ 119 120 typedef struct { 121 PetscReal gravity; 122 } ShallowCtx; 123 124 static inline void ShallowFlux(ShallowCtx *phys,const PetscScalar *u,PetscScalar *f) 125 { 126 f[0] = u[1]; 127 f[1] = PetscSqr(u[1])/u[0] + 0.5*phys->gravity*PetscSqr(u[0]); 128 } 129 130 static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 131 { 132 ShallowCtx *phys = (ShallowCtx*)vctx; 133 PetscScalar g = phys->gravity,ustar[2],cL,cR,c,cstar; 134 struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]},star; 135 PetscInt i; 136 137 PetscFunctionBeginUser; 138 PetscCheck((L.h > 0 && R.h > 0),PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); 139 cL = PetscSqrtScalar(g*L.h); 140 cR = PetscSqrtScalar(g*R.h); 141 c = PetscMax(cL,cR); 142 { 143 /* Solve for star state */ 144 const PetscInt maxits = 50; 145 PetscScalar tmp,res,res0=0,h0,h = 0.5*(L.h + R.h); /* initial guess */ 146 h0 = h; 147 for (i=0; i<maxits; i++) { 148 PetscScalar fr,fl,dfr,dfl; 149 fl = (L.h < h) 150 ? PetscSqrtScalar(0.5*g*(h*h - L.h*L.h)*(1/L.h - 1/h)) /* shock */ 151 : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*L.h); /* rarefaction */ 152 fr = (R.h < h) 153 ? PetscSqrtScalar(0.5*g*(h*h - R.h*R.h)*(1/R.h - 1/h)) /* shock */ 154 : 2*PetscSqrtScalar(g*h) - 2*PetscSqrtScalar(g*R.h); /* rarefaction */ 155 res = R.u - L.u + fr + fl; 156 PetscCheck(!PetscIsInfOrNanScalar(res),PETSC_COMM_SELF,PETSC_ERR_FP,"Infinity or Not-a-Number generated in computation"); 157 if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h-h0) < PETSC_SQRT_MACHINE_EPSILON)) { 158 star.h = h; 159 star.u = L.u - fl; 160 goto converged; 161 } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */ 162 h = 0.8*h0 + 0.2*h; 163 continue; 164 } 165 /* Accept the last step and take another */ 166 res0 = res; 167 h0 = h; 168 dfl = (L.h < h) ? 0.5/fl*0.5*g*(-L.h*L.h/(h*h) - 1 + 2*h/L.h) : PetscSqrtScalar(g/h); 169 dfr = (R.h < h) ? 0.5/fr*0.5*g*(-R.h*R.h/(h*h) - 1 + 2*h/R.h) : PetscSqrtScalar(g/h); 170 tmp = h - res/(dfr+dfl); 171 if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */ 172 else h = tmp; 173 PetscCheck(((h > 0) && PetscIsNormalScalar(h)),PETSC_COMM_SELF,PETSC_ERR_FP,"non-normal iterate h=%g",(double)h); 174 } 175 SETERRQ(PETSC_COMM_SELF,1,"Newton iteration for star.h diverged after %" PetscInt_FMT " iterations",i); 176 } 177 converged: 178 cstar = PetscSqrtScalar(g*star.h); 179 if (L.u-cL < 0 && 0 < star.u-cstar) { /* 1-wave is sonic rarefaction */ 180 PetscScalar ufan[2]; 181 ufan[0] = 1/g*PetscSqr(L.u/3 + 2./3*cL); 182 ufan[1] = PetscSqrtScalar(g*ufan[0])*ufan[0]; 183 ShallowFlux(phys,ufan,flux); 184 } else if (star.u+cstar < 0 && 0 < R.u+cR) { /* 2-wave is sonic rarefaction */ 185 PetscScalar ufan[2]; 186 ufan[0] = 1/g*PetscSqr(R.u/3 - 2./3*cR); 187 ufan[1] = -PetscSqrtScalar(g*ufan[0])*ufan[0]; 188 ShallowFlux(phys,ufan,flux); 189 } else if ((L.h >= star.h && L.u-c >= 0) || (L.h<star.h && (star.h*star.u-L.h*L.u)/(star.h-L.h) > 0)) { 190 /* 1-wave is right-travelling shock (supersonic) */ 191 ShallowFlux(phys,uL,flux); 192 } else if ((star.h <= R.h && R.u+c <= 0) || (star.h>R.h && (R.h*R.u-star.h*star.h)/(R.h-star.h) < 0)) { 193 /* 2-wave is left-travelling shock (supersonic) */ 194 ShallowFlux(phys,uR,flux); 195 } else { 196 ustar[0] = star.h; 197 ustar[1] = star.h*star.u; 198 ShallowFlux(phys,ustar,flux); 199 } 200 *maxspeed = MaxAbs(MaxAbs(star.u-cstar,star.u+cstar),MaxAbs(L.u-cL,R.u+cR)); 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 205 { 206 ShallowCtx *phys = (ShallowCtx*)vctx; 207 PetscScalar g = phys->gravity,fL[2],fR[2],s; 208 struct {PetscScalar h,u;} L = {uL[0],uL[1]/uL[0]},R = {uR[0],uR[1]/uR[0]}; 209 PetscReal tol = 1e-6; 210 211 PetscFunctionBeginUser; 212 /* Positivity preserving modification*/ 213 if (L.h < tol) L.u = 0.0; 214 if (R.h < tol) R.u = 0.0; 215 216 /*simple pos preserve limiter*/ 217 if (L.h < 0) L.h = 0; 218 if (R.h < 0) R.h = 0; 219 220 ShallowFlux(phys,uL,fL); 221 ShallowFlux(phys,uR,fR); 222 223 s = PetscMax(PetscAbs(L.u)+PetscSqrtScalar(g*L.h),PetscAbs(R.u)+PetscSqrtScalar(g*R.h)); 224 flux[0] = 0.5*(fL[0] + fR[0]) + 0.5*s*(L.h - R.h); 225 flux[1] = 0.5*(fL[1] + fR[1]) + 0.5*s*(uL[1] - uR[1]); 226 *maxspeed = s; 227 PetscFunctionReturn(0); 228 } 229 230 static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 231 { 232 PetscInt i,j; 233 234 PetscFunctionBeginUser; 235 for (i=0; i<m; i++) { 236 for (j=0; j<m; j++) Xi[i*m+j] = X[i*m+j] = (PetscScalar)(i==j); 237 speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */ 238 } 239 PetscFunctionReturn(0); 240 } 241 242 static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 243 { 244 ShallowCtx *phys = (ShallowCtx*)vctx; 245 PetscReal c; 246 PetscReal tol = 1e-6; 247 248 PetscFunctionBeginUser; 249 c = PetscSqrtScalar(u[0]*phys->gravity); 250 251 if (u[0] < tol) { /*Use conservative variables*/ 252 X[0*2+0] = 1; 253 X[0*2+1] = 0; 254 X[1*2+0] = 0; 255 X[1*2+1] = 1; 256 } else { 257 speeds[0] = u[1]/u[0] - c; 258 speeds[1] = u[1]/u[0] + c; 259 X[0*2+0] = 1; 260 X[0*2+1] = speeds[0]; 261 X[1*2+0] = 1; 262 X[1*2+1] = speeds[1]; 263 } 264 265 PetscCall(PetscArraycpy(Xi,X,4)); 266 PetscCall(PetscKernel_A_gets_inverse_A_2(Xi,0,PETSC_FALSE,NULL)); 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode PhysicsSample_Shallow(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 271 { 272 PetscFunctionBeginUser; 273 PetscCheck(t <= 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Exact solutions not implemented for t > 0"); 274 switch (initial) { 275 case 0: 276 u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */ 277 u[1] = (x < 0) ? 0 : 0; 278 break; 279 case 1: 280 u[0] = (x < 10) ? 1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */ 281 u[1] = (x < 10) ? 2.5 : 0; 282 break; 283 case 2: 284 u[0] = (x < 25) ? 1 : 1; 285 u[1] = (x < 25) ? -5 : 5; 286 break; 287 case 3: 288 u[0] = (x < 20) ? 1 : 1e-12; 289 u[1] = (x < 20) ? 0 : 0; 290 break; 291 case 4: 292 u[0] = (x < 30) ? 1e-12 : 1; 293 u[1] = (x < 30) ? 0 : 0; 294 break; 295 case 5: 296 u[0] = (x < 25) ? 0.1 : 0.1; 297 u[1] = (x < 25) ? -0.3 : 0.3; 298 break; 299 case 6: 300 u[0] = 1+0.5*PetscSinReal(2*PETSC_PI*x); 301 u[1] = 1*u[0]; 302 break; 303 case 7: 304 if (x < -0.1) { 305 u[0] = 1e-9; 306 u[1] = 0.0; 307 } else if (x < 0.1) { 308 u[0] = 1.0; 309 u[1] = 0.0; 310 } else { 311 u[0] = 1e-9; 312 u[1] = 0.0; 313 } 314 break; 315 case 8: 316 if (x < -0.1) { 317 u[0] = 2; 318 u[1] = 0.0; 319 } else if (x < 0.1) { 320 u[0] = 3.0; 321 u[1] = 0.0; 322 } else { 323 u[0] = 2; 324 u[1] = 0.0; 325 } 326 break; 327 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 328 } 329 PetscFunctionReturn(0); 330 } 331 332 /* Implements inflow conditions for the given initial conditions. Which conditions are actually enforced depends on 333 on the results of PhysicsSetInflowType_Shallow. */ 334 static PetscErrorCode PhysicsInflow_Shallow(void *vctx,PetscReal t,PetscReal x,PetscReal *u) 335 { 336 FVCtx *ctx = (FVCtx*)vctx; 337 338 PetscFunctionBeginUser; 339 if (ctx->bctype == FVBC_INFLOW) { 340 switch (ctx->initial) { 341 case 0: 342 case 1: 343 case 2: 344 case 3: 345 case 4: 346 case 5: 347 u[0] = 0; u[1] = 0.0; /* Left boundary conditions */ 348 u[2] = 0; u[3] = 0.0; /* Right boundary conditions */ 349 break; 350 case 6: 351 u[0] = 0; u[1] = 0.0; /* Left boundary conditions */ 352 u[2] = 0; u[3] = 0.0; /* Right boundary conditions */ 353 break; 354 case 7: 355 u[0] = 0; u[1] = 0.0; /* Left boundary conditions */ 356 u[2] = 0; u[3] = 0.0; /* Right boundary conditions */ 357 break; 358 case 8: 359 u[0] = 0; u[1] = 1.0; /* Left boundary conditions */ 360 u[2] = 0; u[3] = -1.0;/* Right boundary conditions */ 361 break; 362 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 363 } 364 } 365 PetscFunctionReturn(0); 366 } 367 368 /* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */ 369 static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx) 370 { 371 PetscFunctionBeginUser; 372 switch (ctx->initial) { 373 case 0: 374 case 1: 375 case 2: 376 case 3: 377 case 4: 378 case 5: 379 case 6: 380 case 7: 381 case 8: /* Fix left and right momentum, height is outflow */ 382 ctx->physics2.bcinflowindex[0] = PETSC_FALSE; 383 ctx->physics2.bcinflowindex[1] = PETSC_TRUE; 384 ctx->physics2.bcinflowindex[2] = PETSC_FALSE; 385 ctx->physics2.bcinflowindex[3] = PETSC_TRUE; 386 break; 387 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 388 } 389 PetscFunctionReturn(0); 390 } 391 392 static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx) 393 { 394 ShallowCtx *user; 395 PetscFunctionList rlist = 0,rclist = 0; 396 char rname[256] = "rusanov",rcname[256] = "characteristic"; 397 398 PetscFunctionBeginUser; 399 PetscCall(PetscNew(&user)); 400 ctx->physics2.sample2 = PhysicsSample_Shallow; 401 ctx->physics2.inflow = PhysicsInflow_Shallow; 402 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; 403 ctx->physics2.riemann2 = PhysicsRiemann_Shallow_Rusanov; 404 ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow; 405 ctx->physics2.user = user; 406 ctx->physics2.dof = 2; 407 408 PetscMalloc1(2*(ctx->physics2.dof),&ctx->physics2.bcinflowindex); 409 PhysicsSetInflowType_Shallow(ctx); 410 411 PetscCall(PetscStrallocpy("height",&ctx->physics2.fieldname[0])); 412 PetscCall(PetscStrallocpy("momentum",&ctx->physics2.fieldname[1])); 413 414 user->gravity = 9.81; 415 416 PetscCall(RiemannListAdd_2WaySplit(&rlist,"exact", PhysicsRiemann_Shallow_Exact)); 417 PetscCall(RiemannListAdd_2WaySplit(&rlist,"rusanov",PhysicsRiemann_Shallow_Rusanov)); 418 PetscCall(ReconstructListAdd_2WaySplit(&rclist,"characteristic",PhysicsCharacteristic_Shallow)); 419 PetscCall(ReconstructListAdd_2WaySplit(&rclist,"conservative",PhysicsCharacteristic_Conservative)); 420 PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for Shallow",""); 421 PetscCall(PetscOptionsReal("-physics_shallow_gravity","Gravity","",user->gravity,&user->gravity,NULL)); 422 PetscCall(PetscOptionsFList("-physics_shallow_riemann","Riemann solver","",rlist,rname,rname,sizeof(rname),NULL)); 423 PetscCall(PetscOptionsFList("-physics_shallow_reconstruct","Reconstruction","",rclist,rcname,rcname,sizeof(rcname),NULL)); 424 PetscOptionsEnd(); 425 PetscCall(RiemannListFind_2WaySplit(rlist,rname,&ctx->physics2.riemann2)); 426 PetscCall(ReconstructListFind_2WaySplit(rclist,rcname,&ctx->physics2.characteristic2)); 427 PetscCall(PetscFunctionListDestroy(&rlist)); 428 PetscCall(PetscFunctionListDestroy(&rclist)); 429 PetscFunctionReturn(0); 430 } 431 432 PetscErrorCode FVSample_2WaySplit(FVCtx *ctx,DM da,PetscReal time,Vec U) 433 { 434 PetscScalar *u,*uj,xj,xi; 435 PetscInt i,j,k,dof,xs,xm,Mx; 436 const PetscInt N = 200; 437 PetscReal hs,hf; 438 439 PetscFunctionBeginUser; 440 PetscCheck(ctx->physics2.sample2,PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 441 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 442 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 443 PetscCall(DMDAVecGetArray(da,U,&u)); 444 PetscCall(PetscMalloc1(dof,&uj)); 445 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 446 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 447 for (i=xs; i<xs+xm; i++) { 448 if (i < ctx->sf) { 449 xi = ctx->xmin+0.5*hs+i*hs; 450 /* Integrate over cell i using trapezoid rule with N points. */ 451 for (k=0; k<dof; k++) u[i*dof+k] = 0; 452 for (j=0; j<N+1; j++) { 453 xj = xi+hs*(j-N/2)/(PetscReal)N; 454 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 455 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 456 } 457 } else if (i < ctx->fs) { 458 xi = ctx->xmin+ctx->sf*hs+0.5*hf+(i-ctx->sf)*hf; 459 /* Integrate over cell i using trapezoid rule with N points. */ 460 for (k=0; k<dof; k++) u[i*dof+k] = 0; 461 for (j=0; j<N+1; j++) { 462 xj = xi+hf*(j-N/2)/(PetscReal)N; 463 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 464 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 465 } 466 } else { 467 xi = ctx->xmin+ctx->sf*hs+(ctx->fs-ctx->sf)*hf+0.5*hs+(i-ctx->fs)*hs; 468 /* Integrate over cell i using trapezoid rule with N points. */ 469 for (k=0; k<dof; k++) u[i*dof+k] = 0; 470 for (j=0; j<N+1; j++) { 471 xj = xi+hs*(j-N/2)/(PetscReal)N; 472 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 473 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 474 } 475 } 476 } 477 PetscCall(DMDAVecRestoreArray(da,U,&u)); 478 PetscCall(PetscFree(uj)); 479 PetscFunctionReturn(0); 480 } 481 482 static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1) 483 { 484 Vec Y; 485 PetscInt i,Mx; 486 const PetscScalar *ptr_X,*ptr_Y; 487 PetscReal hs,hf; 488 489 PetscFunctionBeginUser; 490 PetscCall(VecGetSize(X,&Mx)); 491 PetscCall(VecDuplicate(X,&Y)); 492 PetscCall(FVSample_2WaySplit(ctx,da,t,Y)); 493 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 494 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 495 PetscCall(VecGetArrayRead(X,&ptr_X)); 496 PetscCall(VecGetArrayRead(Y,&ptr_Y)); 497 for (i=0; i<Mx; i++) { 498 if (i < ctx->sf || i > ctx->fs -1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]); 499 else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]); 500 } 501 PetscCall(VecRestoreArrayRead(X,&ptr_X)); 502 PetscCall(VecRestoreArrayRead(Y,&ptr_Y)); 503 PetscCall(VecDestroy(&Y)); 504 PetscFunctionReturn(0); 505 } 506 507 PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 508 { 509 FVCtx *ctx = (FVCtx*)vctx; 510 PetscInt i,j,k,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs; 511 PetscReal hxf,hxs,cfl_idt = 0; 512 PetscScalar *x,*f,*slope; 513 Vec Xloc; 514 DM da; 515 516 PetscFunctionBeginUser; 517 PetscCall(TSGetDM(ts,&da)); 518 PetscCall(DMGetLocalVector(da,&Xloc)); /* Xloc contains ghost points */ 519 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); /* Mx is the number of center points */ 520 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 521 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 522 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); /* X is solution vector which does not contain ghost points */ 523 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 524 525 PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */ 526 527 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 528 PetscCall(DMDAVecGetArray(da,F,&f)); 529 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); /* contains ghost points */ 530 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 531 532 if (ctx->bctype == FVBC_OUTFLOW) { 533 for (i=xs-2; i<0; i++) { 534 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 535 } 536 for (i=Mx; i<xs+xm+2; i++) { 537 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 538 } 539 } 540 541 if (ctx->bctype == FVBC_INFLOW) { 542 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 543 pages 137-138 for the scheme. */ 544 if (xs==0) { /* Left Boundary */ 545 ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub); 546 for (j=0; j<dof; j++) { 547 if (ctx->physics2.bcinflowindex[j]) { 548 for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j]; 549 } else { 550 for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */ 551 } 552 } 553 } 554 if (xs+xm==Mx) { /* Right Boundary */ 555 ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub); 556 for (j=0; j<dof; j++) { 557 if (ctx->physics2.bcinflowindex[dof+j]) { 558 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j]; 559 } else { 560 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */ 561 } 562 } 563 } 564 } 565 566 for (i=xs-1; i<xs+xm+1; i++) { 567 struct _LimitInfo info; 568 PetscScalar *cjmpL,*cjmpR; 569 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 570 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 571 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 572 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 573 cjmpL = &ctx->cjmpLR[0]; 574 cjmpR = &ctx->cjmpLR[dof]; 575 for (j=0; j<dof; j++) { 576 PetscScalar jmpL,jmpR; 577 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 578 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 579 for (k=0; k<dof; k++) { 580 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 581 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 582 } 583 } 584 /* Apply limiter to the left and right characteristic jumps */ 585 info.m = dof; 586 info.hxs = hxs; 587 info.hxf = hxf; 588 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 589 for (j=0; j<dof; j++) { 590 PetscScalar tmp = 0; 591 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 592 slope[i*dof+j] = tmp; 593 } 594 } 595 596 for (i=xs; i<xs+xm+1; i++) { 597 PetscReal maxspeed; 598 PetscScalar *uL,*uR; 599 uL = &ctx->uLR[0]; 600 uR = &ctx->uLR[dof]; 601 if (i < sf) { /* slow region */ 602 for (j=0; j<dof; j++) { 603 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 604 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 605 } 606 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 607 if (i > xs) { 608 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 609 } 610 if (i < xs+xm) { 611 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 612 } 613 } else if (i == sf) { /* interface between the slow region and the fast region */ 614 for (j=0; j<dof; j++) { 615 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 616 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 617 } 618 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 619 if (i > xs) { 620 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 621 } 622 if (i < xs+xm) { 623 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf; 624 } 625 } else if (i > sf && i < fs) { /* fast region */ 626 for (j=0; j<dof; j++) { 627 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 628 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 629 } 630 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 631 if (i > xs) { 632 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf; 633 } 634 if (i < xs+xm) { 635 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf; 636 } 637 } else if (i == fs) { /* interface between the fast region and the slow region */ 638 for (j=0; j<dof; j++) { 639 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 640 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 641 } 642 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 643 if (i > xs) { 644 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf; 645 } 646 if (i < xs+xm) { 647 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 648 } 649 } else { /* slow region */ 650 for (j=0; j<dof; j++) { 651 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 652 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 653 } 654 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 655 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */ 656 if (i > xs) { 657 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 658 } 659 if (i < xs+xm) { 660 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 661 } 662 } 663 } 664 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 665 PetscCall(DMDAVecRestoreArray(da,F,&f)); 666 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 667 PetscCall(DMRestoreLocalVector(da,&Xloc)); 668 PetscCallMPI(MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da))); 669 if (0) { 670 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 671 PetscReal dt,tnow; 672 PetscCall(TSGetTimeStep(ts,&dt)); 673 PetscCall(TSGetTime(ts,&tnow)); 674 if (dt > 0.5/ctx->cfl_idt) { 675 if (1) { 676 PetscCall(PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt))); 677 } else SETERRQ(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt)); 678 } 679 } 680 PetscFunctionReturn(0); 681 } 682 683 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */ 684 PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 685 { 686 FVCtx *ctx = (FVCtx*)vctx; 687 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth; 688 PetscReal hxs,hxf,cfl_idt = 0; 689 PetscScalar *x,*f,*slope; 690 Vec Xloc; 691 DM da; 692 693 PetscFunctionBeginUser; 694 PetscCall(TSGetDM(ts,&da)); 695 PetscCall(DMGetLocalVector(da,&Xloc)); 696 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 697 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 698 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 699 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 700 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 701 PetscCall(VecZeroEntries(F)); 702 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 703 PetscCall(VecGetArray(F,&f)); 704 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 705 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 706 707 if (ctx->bctype == FVBC_OUTFLOW) { 708 for (i=xs-2; i<0; i++) { 709 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 710 } 711 for (i=Mx; i<xs+xm+2; i++) { 712 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 713 } 714 } 715 if (ctx->bctype == FVBC_INFLOW) { 716 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 717 pages 137-138 for the scheme. */ 718 if (xs==0) { /* Left Boundary */ 719 ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub); 720 for (j=0; j<dof; j++) { 721 if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) { 722 for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j]; 723 } else { 724 for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */ 725 } 726 } 727 } 728 if (xs+xm==Mx) { /* Right Boundary */ 729 ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub); 730 for (j=0; j<dof; j++) { 731 if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) { 732 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j]; 733 } else { 734 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */ 735 } 736 } 737 } 738 } 739 for (i=xs-1; i<xs+xm+1; i++) { 740 struct _LimitInfo info; 741 PetscScalar *cjmpL,*cjmpR; 742 if (i < sf-lsbwidth+1 || i > fs+rsbwidth-2) { /* slow components and the first and last fast components */ 743 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 744 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 745 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 746 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 747 cjmpL = &ctx->cjmpLR[0]; 748 cjmpR = &ctx->cjmpLR[dof]; 749 for (j=0; j<dof; j++) { 750 PetscScalar jmpL,jmpR; 751 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 752 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 753 for (k=0; k<dof; k++) { 754 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 755 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 756 } 757 } 758 /* Apply limiter to the left and right characteristic jumps */ 759 info.m = dof; 760 info.hxs = hxs; 761 info.hxf = hxf; 762 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 763 for (j=0; j<dof; j++) { 764 PetscScalar tmp = 0; 765 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 766 slope[i*dof+j] = tmp; 767 } 768 } 769 } 770 771 for (i=xs; i<xs+xm+1; i++) { 772 PetscReal maxspeed; 773 PetscScalar *uL,*uR; 774 uL = &ctx->uLR[0]; 775 uR = &ctx->uLR[dof]; 776 if (i < sf-lsbwidth) { /* slow region */ 777 for (j=0; j<dof; j++) { 778 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 779 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 780 } 781 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 782 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */ 783 if (i > xs) { 784 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 785 } 786 if (i < xs+xm) { 787 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 788 islow++; 789 } 790 } 791 if (i == sf-lsbwidth) { /* interface between the slow region and the fast region */ 792 for (j=0; j<dof; j++) { 793 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 794 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 795 } 796 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 797 if (i > xs) { 798 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 799 } 800 } 801 if (i == fs+rsbwidth) { /* slow region */ 802 for (j=0; j<dof; j++) { 803 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 804 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 805 } 806 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 807 if (i < xs+xm) { 808 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 809 islow++; 810 } 811 } 812 if (i > fs+rsbwidth) { /* slow region */ 813 for (j=0; j<dof; j++) { 814 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 815 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 816 } 817 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 818 if (i > xs) { 819 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 820 } 821 if (i < xs+xm) { 822 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 823 islow++; 824 } 825 } 826 } 827 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 828 PetscCall(VecRestoreArray(F,&f)); 829 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 830 PetscCall(DMRestoreLocalVector(da,&Xloc)); 831 PetscCallMPI(MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da))); 832 PetscFunctionReturn(0); 833 } 834 835 PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 836 { 837 FVCtx *ctx = (FVCtx*)vctx; 838 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth; 839 PetscReal hxs,hxf; 840 PetscScalar *x,*f,*slope; 841 Vec Xloc; 842 DM da; 843 844 PetscFunctionBeginUser; 845 PetscCall(TSGetDM(ts,&da)); 846 PetscCall(DMGetLocalVector(da,&Xloc)); 847 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 848 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 849 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 850 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 851 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 852 PetscCall(VecZeroEntries(F)); 853 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 854 PetscCall(VecGetArray(F,&f)); 855 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 856 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 857 858 if (ctx->bctype == FVBC_OUTFLOW) { 859 for (i=xs-2; i<0; i++) { 860 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 861 } 862 for (i=Mx; i<xs+xm+2; i++) { 863 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 864 } 865 } 866 if (ctx->bctype == FVBC_INFLOW) { 867 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 868 pages 137-138 for the scheme. */ 869 if (xs==0) { /* Left Boundary */ 870 ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub); 871 for (j=0; j<dof; j++) { 872 if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) { 873 for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j]; 874 } else { 875 for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */ 876 } 877 } 878 } 879 if (xs+xm==Mx) { /* Right Boundary */ 880 ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub); 881 for (j=0; j<dof; j++) { 882 if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) { 883 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j]; 884 } else { 885 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */ 886 } 887 } 888 } 889 } 890 for (i=xs-1; i<xs+xm+1; i++) { 891 struct _LimitInfo info; 892 PetscScalar *cjmpL,*cjmpR; 893 if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) { 894 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 895 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 896 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 897 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 898 cjmpL = &ctx->cjmpLR[0]; 899 cjmpR = &ctx->cjmpLR[dof]; 900 for (j=0; j<dof; j++) { 901 PetscScalar jmpL,jmpR; 902 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 903 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 904 for (k=0; k<dof; k++) { 905 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 906 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 907 } 908 } 909 /* Apply limiter to the left and right characteristic jumps */ 910 info.m = dof; 911 info.hxs = hxs; 912 info.hxf = hxf; 913 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 914 for (j=0; j<dof; j++) { 915 PetscScalar tmp = 0; 916 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 917 slope[i*dof+j] = tmp; 918 } 919 } 920 } 921 922 for (i=xs; i<xs+xm+1; i++) { 923 PetscReal maxspeed; 924 PetscScalar *uL,*uR; 925 uL = &ctx->uLR[0]; 926 uR = &ctx->uLR[dof]; 927 if (i == sf-lsbwidth) { 928 for (j=0; j<dof; j++) { 929 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 930 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 931 } 932 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 933 if (i < xs+xm) { 934 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 935 islow++; 936 } 937 } 938 if (i > sf-lsbwidth && i < sf) { 939 for (j=0; j<dof; j++) { 940 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 941 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 942 } 943 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 944 if (i > xs) { 945 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 946 } 947 if (i < xs+xm) { 948 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 949 islow++; 950 } 951 } 952 if (i == sf) { /* interface between the slow region and the fast region */ 953 for (j=0; j<dof; j++) { 954 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 955 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 956 } 957 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 958 if (i > xs) { 959 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 960 } 961 } 962 if (i == fs) { /* interface between the fast region and the slow region */ 963 for (j=0; j<dof; j++) { 964 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 965 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 966 } 967 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 968 if (i < xs+xm) { 969 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 970 islow++; 971 } 972 } 973 if (i > fs && i < fs+rsbwidth) { 974 for (j=0; j<dof; j++) { 975 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 976 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 977 } 978 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 979 if (i > xs) { 980 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 981 } 982 if (i < xs+xm) { 983 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 984 islow++; 985 } 986 } 987 if (i == fs+rsbwidth) { 988 for (j=0; j<dof; j++) { 989 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 990 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 991 } 992 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 993 if (i > xs) { 994 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 995 } 996 } 997 } 998 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 999 PetscCall(VecRestoreArray(F,&f)); 1000 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 1001 PetscCall(DMRestoreLocalVector(da,&Xloc)); 1002 PetscFunctionReturn(0); 1003 } 1004 1005 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */ 1006 PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 1007 { 1008 FVCtx *ctx = (FVCtx*)vctx; 1009 PetscInt i,j,k,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs; 1010 PetscReal hxs,hxf; 1011 PetscScalar *x,*f,*slope; 1012 Vec Xloc; 1013 DM da; 1014 1015 PetscFunctionBeginUser; 1016 PetscCall(TSGetDM(ts,&da)); 1017 PetscCall(DMGetLocalVector(da,&Xloc)); 1018 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 1019 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 1020 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 1021 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 1022 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 1023 PetscCall(VecZeroEntries(F)); 1024 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 1025 PetscCall(VecGetArray(F,&f)); 1026 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 1027 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 1028 1029 if (ctx->bctype == FVBC_OUTFLOW) { 1030 for (i=xs-2; i<0; i++) { 1031 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 1032 } 1033 for (i=Mx; i<xs+xm+2; i++) { 1034 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 1035 } 1036 } 1037 if (ctx->bctype == FVBC_INFLOW) { 1038 /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253 1039 pages 137-138 for the scheme. */ 1040 if (xs==0) { /* Left Boundary */ 1041 ctx->physics2.inflow(ctx,time,ctx->xmin,ctx->ub); 1042 for (j=0; j<dof; j++) { 1043 if (ctx->physics2.bcinflowindex[j]==PETSC_TRUE) { 1044 for (i=-2; i<0; i++) x[i*dof+j] = 2.0*ctx->ub[j]-x[-(i+1)*dof+j]; 1045 } else { 1046 for (i=-2; i<0; i++) x[i*dof+j] = x[j]; /* Outflow */ 1047 } 1048 } 1049 } 1050 if (xs+xm==Mx) { /* Right Boundary */ 1051 ctx->physics2.inflow(ctx,time,ctx->xmax,ctx->ub); 1052 for (j=0; j<dof; j++) { 1053 if (ctx->physics2.bcinflowindex[dof+j]==PETSC_TRUE) { 1054 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = 2.0*ctx->ub[dof+j]-x[(2*Mx-(i+1))*dof+j]; 1055 } else { 1056 for (i=Mx; i<Mx+2; i++) x[i*dof+j] = x[(Mx-1)*dof+j]; /* Outflow */ 1057 } 1058 } 1059 } 1060 } 1061 for (i=xs-1; i<xs+xm+1; i++) { 1062 struct _LimitInfo info; 1063 PetscScalar *cjmpL,*cjmpR; 1064 if (i > sf-2 && i < fs+1) { 1065 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 1066 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 1067 cjmpL = &ctx->cjmpLR[0]; 1068 cjmpR = &ctx->cjmpLR[dof]; 1069 for (j=0; j<dof; j++) { 1070 PetscScalar jmpL,jmpR; 1071 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 1072 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 1073 for (k=0; k<dof; k++) { 1074 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 1075 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 1076 } 1077 } 1078 /* Apply limiter to the left and right characteristic jumps */ 1079 info.m = dof; 1080 info.hxs = hxs; 1081 info.hxf = hxf; 1082 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 1083 for (j=0; j<dof; j++) { 1084 PetscScalar tmp = 0; 1085 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 1086 slope[i*dof+j] = tmp; 1087 } 1088 } 1089 } 1090 1091 for (i=xs; i<xs+xm+1; i++) { 1092 PetscReal maxspeed; 1093 PetscScalar *uL,*uR; 1094 uL = &ctx->uLR[0]; 1095 uR = &ctx->uLR[dof]; 1096 if (i == sf) { /* interface between the slow region and the fast region */ 1097 for (j=0; j<dof; j++) { 1098 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 1099 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 1100 } 1101 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 1102 if (i < xs+xm) { 1103 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf; 1104 ifast++; 1105 } 1106 } 1107 if (i > sf && i < fs) { /* fast region */ 1108 for (j=0; j<dof; j++) { 1109 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 1110 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 1111 } 1112 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 1113 if (i > xs) { 1114 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf; 1115 } 1116 if (i < xs+xm) { 1117 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf; 1118 ifast++; 1119 } 1120 } 1121 if (i == fs) { /* interface between the fast region and the slow region */ 1122 for (j=0; j<dof; j++) { 1123 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 1124 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 1125 } 1126 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 1127 if (i > xs) { 1128 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf; 1129 } 1130 } 1131 } 1132 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 1133 PetscCall(VecRestoreArray(F,&f)); 1134 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 1135 PetscCall(DMRestoreLocalVector(da,&Xloc)); 1136 PetscFunctionReturn(0); 1137 } 1138 1139 int main(int argc,char *argv[]) 1140 { 1141 char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m"; 1142 PetscFunctionList limiters = 0,physics = 0; 1143 MPI_Comm comm; 1144 TS ts; 1145 DM da; 1146 Vec X,X0,R; 1147 FVCtx ctx; 1148 PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast =0,islowbuffer = 0,*index_slow,*index_fast,*index_slowbuffer; 1149 PetscBool view_final = PETSC_FALSE; 1150 PetscReal ptime,maxtime; 1151 1152 PetscCall(PetscInitialize(&argc,&argv,0,help)); 1153 comm = PETSC_COMM_WORLD; 1154 PetscCall(PetscMemzero(&ctx,sizeof(ctx))); 1155 1156 /* Register limiters to be available on the command line */ 1157 PetscCall(PetscFunctionListAdd(&limiters,"upwind" ,Limit2_Upwind)); 1158 PetscCall(PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit2_LaxWendroff)); 1159 PetscCall(PetscFunctionListAdd(&limiters,"beam-warming" ,Limit2_BeamWarming)); 1160 PetscCall(PetscFunctionListAdd(&limiters,"fromm" ,Limit2_Fromm)); 1161 PetscCall(PetscFunctionListAdd(&limiters,"minmod" ,Limit2_Minmod)); 1162 PetscCall(PetscFunctionListAdd(&limiters,"superbee" ,Limit2_Superbee)); 1163 PetscCall(PetscFunctionListAdd(&limiters,"mc" ,Limit2_MC)); 1164 PetscCall(PetscFunctionListAdd(&limiters,"koren3" ,Limit2_Koren3)); 1165 1166 /* Register physical models to be available on the command line */ 1167 PetscCall(PetscFunctionListAdd(&physics,"shallow" ,PhysicsCreate_Shallow)); 1168 PetscCall(PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect)); 1169 1170 ctx.comm = comm; 1171 ctx.cfl = 0.9; 1172 ctx.bctype = FVBC_PERIODIC; 1173 ctx.xmin = -1.0; 1174 ctx.xmax = 1.0; 1175 ctx.initial = 1; 1176 ctx.hratio = 2; 1177 maxtime = 10.0; 1178 ctx.simulation = PETSC_FALSE; 1179 PetscOptionsBegin(comm,NULL,"Finite Volume solver options",""); 1180 PetscCall(PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL)); 1181 PetscCall(PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL)); 1182 PetscCall(PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL)); 1183 PetscCall(PetscOptionsFList("-physics","Name of physics model to use","",physics,physname,physname,sizeof(physname),NULL)); 1184 PetscCall(PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL)); 1185 PetscCall(PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final)); 1186 PetscCall(PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL)); 1187 PetscCall(PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL)); 1188 PetscCall(PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL)); 1189 PetscCall(PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL)); 1190 PetscCall(PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL)); 1191 PetscCall(PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL)); 1192 PetscOptionsEnd(); 1193 1194 /* Choose the limiter from the list of registered limiters */ 1195 PetscCall(PetscFunctionListFind(limiters,lname,&ctx.limit2)); 1196 PetscCheck(ctx.limit2,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Limiter '%s' not found",lname); 1197 1198 /* Choose the physics from the list of registered models */ 1199 { 1200 PetscErrorCode (*r)(FVCtx*); 1201 PetscCall(PetscFunctionListFind(physics,physname,&r)); 1202 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 1203 /* Create the physics, will set the number of fields and their names */ 1204 PetscCall((*r)(&ctx)); 1205 } 1206 1207 /* Create a DMDA to manage the parallel grid */ 1208 PetscCall(DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da)); 1209 PetscCall(DMSetFromOptions(da)); 1210 PetscCall(DMSetUp(da)); 1211 /* Inform the DMDA of the field names provided by the physics. */ 1212 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 1213 for (i=0; i<ctx.physics2.dof; i++) { 1214 PetscCall(DMDASetFieldName(da,i,ctx.physics2.fieldname[i])); 1215 } 1216 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 1217 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 1218 1219 /* Set coordinates of cell centers */ 1220 PetscCall(DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0)); 1221 1222 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 1223 PetscCall(PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope)); 1224 PetscCall(PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds)); 1225 PetscCall(PetscMalloc1(2*dof,&ctx.ub)); 1226 1227 /* Create a vector to store the solution and to save the initial state */ 1228 PetscCall(DMCreateGlobalVector(da,&X)); 1229 PetscCall(VecDuplicate(X,&X0)); 1230 PetscCall(VecDuplicate(X,&R)); 1231 1232 /* create index for slow parts and fast parts, 1233 count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */ 1234 count_slow = Mx/(1.0+ctx.hratio/3.0); 1235 PetscCheck(count_slow%2 == 0,PETSC_COMM_WORLD,PETSC_ERR_USER,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hratio/3) is even"); 1236 count_fast = Mx-count_slow; 1237 ctx.sf = count_slow/2; 1238 ctx.fs = ctx.sf+count_fast; 1239 PetscCall(PetscMalloc1(xm*dof,&index_slow)); 1240 PetscCall(PetscMalloc1(xm*dof,&index_fast)); 1241 PetscCall(PetscMalloc1(8*dof,&index_slowbuffer)); 1242 ctx.lsbwidth = 4; 1243 ctx.rsbwidth = 4; 1244 1245 for (i=xs; i<xs+xm; i++) { 1246 if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1) 1247 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 1248 else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1)) 1249 for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k; 1250 else 1251 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 1252 } 1253 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss)); 1254 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf)); 1255 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb)); 1256 1257 /* Create a time-stepping object */ 1258 PetscCall(TSCreate(comm,&ts)); 1259 PetscCall(TSSetDM(ts,da)); 1260 PetscCall(TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx)); 1261 PetscCall(TSRHSSplitSetIS(ts,"slow",ctx.iss)); 1262 PetscCall(TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb)); 1263 PetscCall(TSRHSSplitSetIS(ts,"fast",ctx.isf)); 1264 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx)); 1265 PetscCall(TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx)); 1266 PetscCall(TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx)); 1267 1268 PetscCall(TSSetType(ts,TSMPRK)); 1269 PetscCall(TSSetMaxTime(ts,maxtime)); 1270 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1271 1272 /* Compute initial conditions and starting time step */ 1273 PetscCall(FVSample_2WaySplit(&ctx,da,0,X0)); 1274 PetscCall(FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx)); /* Initial function evaluation, only used to determine max speed */ 1275 PetscCall(VecCopy(X0,X)); /* The function value was not used so we set X=X0 again */ 1276 PetscCall(TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt)); 1277 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 1278 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 1279 { 1280 PetscInt steps; 1281 PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg; 1282 const PetscScalar *ptr_X,*ptr_X0; 1283 const PetscReal hs = (ctx.xmax-ctx.xmin)*3.0/4.0/count_slow; 1284 const PetscReal hf = (ctx.xmax-ctx.xmin)/4.0/count_fast; 1285 1286 PetscCall(TSSolve(ts,X)); 1287 PetscCall(TSGetSolveTime(ts,&ptime)); 1288 PetscCall(TSGetStepNumber(ts,&steps)); 1289 /* calculate the total mass at initial time and final time */ 1290 mass_initial = 0.0; 1291 mass_final = 0.0; 1292 PetscCall(DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0)); 1293 PetscCall(DMDAVecGetArrayRead(da,X,(void*)&ptr_X)); 1294 for (i=xs;i<xs+xm;i++) { 1295 if (i < ctx.sf || i > ctx.fs-1) { 1296 for (k=0; k<dof; k++) { 1297 mass_initial = mass_initial + hs*ptr_X0[i*dof+k]; 1298 mass_final = mass_final + hs*ptr_X[i*dof+k]; 1299 } 1300 } else { 1301 for (k=0; k<dof; k++) { 1302 mass_initial = mass_initial + hf*ptr_X0[i*dof+k]; 1303 mass_final = mass_final + hf*ptr_X[i*dof+k]; 1304 } 1305 } 1306 } 1307 PetscCall(DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0)); 1308 PetscCall(DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X)); 1309 mass_difference = mass_final - mass_initial; 1310 PetscCallMPI(MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm)); 1311 PetscCall(PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg)); 1312 PetscCall(PetscPrintf(comm,"Final time %g, steps %" PetscInt_FMT "\n",(double)ptime,steps)); 1313 PetscCall(PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt))); 1314 if (ctx.exact) { 1315 PetscReal nrm1=0; 1316 PetscCall(SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1)); 1317 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1)); 1318 } 1319 if (ctx.simulation) { 1320 PetscReal nrm1=0; 1321 PetscViewer fd; 1322 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 1323 Vec XR; 1324 PetscBool flg; 1325 const PetscScalar *ptr_XR; 1326 PetscCall(PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg)); 1327 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 1328 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd)); 1329 PetscCall(VecDuplicate(X0,&XR)); 1330 PetscCall(VecLoad(XR,fd)); 1331 PetscCall(PetscViewerDestroy(&fd)); 1332 PetscCall(VecGetArrayRead(X,&ptr_X)); 1333 PetscCall(VecGetArrayRead(XR,&ptr_XR)); 1334 for (i=xs;i<xs+xm;i++) { 1335 if (i < ctx.sf || i > ctx.fs-1) 1336 for (k=0; k<dof; k++) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]); 1337 else 1338 for (k=0; k<dof; k++) nrm1 = nrm1 + hf*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]); 1339 } 1340 PetscCall(VecRestoreArrayRead(X,&ptr_X)); 1341 PetscCall(VecRestoreArrayRead(XR,&ptr_XR)); 1342 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1)); 1343 PetscCall(VecDestroy(&XR)); 1344 } 1345 } 1346 1347 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 1348 if (draw & 0x1) PetscCall(VecView(X0,PETSC_VIEWER_DRAW_WORLD)); 1349 if (draw & 0x2) PetscCall(VecView(X,PETSC_VIEWER_DRAW_WORLD)); 1350 if (draw & 0x4) { 1351 Vec Y; 1352 PetscCall(VecDuplicate(X,&Y)); 1353 PetscCall(FVSample_2WaySplit(&ctx,da,ptime,Y)); 1354 PetscCall(VecAYPX(Y,-1,X)); 1355 PetscCall(VecView(Y,PETSC_VIEWER_DRAW_WORLD)); 1356 PetscCall(VecDestroy(&Y)); 1357 } 1358 1359 if (view_final) { 1360 PetscViewer viewer; 1361 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer)); 1362 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); 1363 PetscCall(VecView(X,viewer)); 1364 PetscCall(PetscViewerPopFormat(viewer)); 1365 PetscCall(PetscViewerDestroy(&viewer)); 1366 } 1367 1368 /* Clean up */ 1369 PetscCall((*ctx.physics2.destroy)(ctx.physics2.user)); 1370 for (i=0; i<ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i])); 1371 PetscCall(PetscFree(ctx.physics2.bcinflowindex)); 1372 PetscCall(PetscFree(ctx.ub)); 1373 PetscCall(PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope)); 1374 PetscCall(PetscFree3(ctx.uLR,ctx.flux,ctx.speeds)); 1375 PetscCall(VecDestroy(&X)); 1376 PetscCall(VecDestroy(&X0)); 1377 PetscCall(VecDestroy(&R)); 1378 PetscCall(DMDestroy(&da)); 1379 PetscCall(TSDestroy(&ts)); 1380 PetscCall(ISDestroy(&ctx.iss)); 1381 PetscCall(ISDestroy(&ctx.isf)); 1382 PetscCall(ISDestroy(&ctx.issb)); 1383 PetscCall(PetscFree(index_slow)); 1384 PetscCall(PetscFree(index_fast)); 1385 PetscCall(PetscFree(index_slowbuffer)); 1386 PetscCall(PetscFunctionListDestroy(&limiters)); 1387 PetscCall(PetscFunctionListDestroy(&physics)); 1388 PetscCall(PetscFinalize()); 1389 return 0; 1390 } 1391 1392 /*TEST 1393 1394 build: 1395 requires: !complex !single 1396 depends: finitevolume1d.c 1397 1398 test: 1399 suffix: 1 1400 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 1401 output_file: output/ex4_1.out 1402 1403 test: 1404 suffix: 2 1405 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 1406 output_file: output/ex4_1.out 1407 1408 test: 1409 suffix: 3 1410 args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0 1411 output_file: output/ex4_3.out 1412 1413 test: 1414 suffix: 4 1415 nsize: 2 1416 args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1 1417 output_file: output/ex4_3.out 1418 1419 test: 1420 suffix: 5 1421 nsize: 4 1422 args: args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_dt 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1 1423 output_file: output/ex4_3.out 1424 TEST*/ 1425