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