Skip to content

Commit

Permalink
fix outflow BCs for the 4th order solver (#2607)
Browse files Browse the repository at this point in the history 
In the interpolation of the interface states for the 4th order solver, we use a one-sided
stencil at the domain boundaries (so long as we are not periodic). This works great for
symmetry, and we further enforce the symmetry by reflecting the interface states.

For outflow BCs, we were trying to be fancy and continue to use this one-sided stencil
but extrapolate into the ghost cell to get a better prediction for the zone just outside of
the boundary. This does not work well and causes a flame problem to do bad things.

This PR removes all the fancy outflow attempts and has us only do a one-sided stencil
for reflecting. I think this is more correct, and will let us deal with inflow eventually,
although we will not be 4th order for outflow. But I am not sure an outflow BC is defined
well enough to understand what it means to be high order there.
  • Loading branch information
@zingale
zingale authored Oct 9, 2023
1 parent 40b9403 commit be6201a
Showing 1 changed file with 18 additions and 84 deletions.
102 changes: 18 additions & 84 deletions Source/hydro/fourth_order.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,52 +37,30 @@ Castro::fourth_interfaces(const Box& bx,
// interpolate to the edges -- this is a_{i-1/2}
// note for non-periodic physical boundaries, we use a special stencil

if (i == domlo[0]+1 && lo_bc[0] != Interior) {
if (i == domlo[0]+1 && lo_bc[0] == Symmetry) {
// use a stencil for the interface that is one zone
// from the left physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i-1,j,k,ncomp) + 13.0_rt*a(i,j,k,ncomp) -
5.0_rt*a(i+1,j,k,ncomp) + a(i+2,j,k,ncomp));

} else if (i == domlo[0] && lo_bc[0] != Interior) {
} else if (i == domlo[0] && lo_bc[0] == Symmetry) {
// use a stencil for when the interface is on the
// left physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i+1,j,k,ncomp) +
13.0_rt*a(i+2,j,k,ncomp) - 3.0_rt*a(i+3,j,k,ncomp));

} else if (i == domhi[0] && hi_bc[0] != Interior) {
} else if (i == domhi[0] && hi_bc[0] == Symmetry) {
// use a stencil for the interface that is one zone
// from the right physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i-1,j,k,ncomp) -
5.0_rt*a(i-2,j,k,ncomp) + a(i-3,j,k,ncomp));

} else if (i == domhi[0]+1 && hi_bc[0] != Interior) {
} else if (i == domhi[0]+1 && hi_bc[0] == Symmetry) {
// use a stencil for when the interface is on the
// right physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i-1,j,k,ncomp) - 23.0_rt*a(i-2,j,k,ncomp) +
13.0_rt*a(i-3,j,k,ncomp) - 3.0_rt*a(i-4,j,k,ncomp));

} else if (i == domlo[0]-1 && lo_bc[0] == Outflow) {
// extrapolate to the domlo[0]-1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(domlo[0],j,k,ncomp) - 6.0_rt*a(domlo[0]+1,j,k,ncomp) + 4.0_rt*a(domlo[0]+2,j,k,ncomp) - a(domlo[0]+3,j,k,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i+1,j,k,ncomp) +
13.0_rt*a(i+2,j,k,ncomp) - 3.0_rt*a(i+3,j,k,ncomp));

} else if (i == domhi[0]+2 && hi_bc[0] == Outflow) {
// extrapolate to the domhi[0]+1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(domhi[0],j,k,ncomp) - 6.0_rt*a(domhi[0]-1,j,k,ncomp) + 4.0_rt*a(domhi[0]-2,j,k,ncomp) - a(domhi[0]-3,j,k,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i-2,j,k,ncomp) +
13.0_rt*a(i-3,j,k,ncomp) - 3.0_rt*a(i-4,j,k,ncomp));

} else {
// regular stencil
a_int(i,j,k) = (7.0_rt/12.0_rt)*(a(i-1,j,k,ncomp) + a(i,j,k,ncomp)) -
Expand All @@ -100,52 +78,30 @@ Castro::fourth_interfaces(const Box& bx,

// interpolate to the edges

if (j == domlo[1]+1 && lo_bc[1] != Interior) {
if (j == domlo[1]+1 && lo_bc[1] == Symmetry) {
// use a stencil for the interface that is one zone
// from the left physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j-1,k,ncomp) + 13.0_rt*a(i,j,k,ncomp) -
5.0_rt*a(i,j+1,k,ncomp) + a(i,j+2,k,ncomp));

} else if (j == domlo[1] && lo_bc[1] != Interior) {
} else if (j == domlo[1] && lo_bc[1] == Symmetry) {
// use a stencil for when the interface is on the
// left physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i,j+1,k,ncomp) +
13.0_rt*a(i,j+2,k,ncomp) - 3.0_rt*a(i,j+3,k,ncomp));

} else if (j == domhi[1] && hi_bc[1] != Interior) {
} else if (j == domhi[1] && hi_bc[1] == Symmetry) {
// use a stencil for the interface that is one zone
// from the right physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i,j-1,k,ncomp) -
5.0_rt*a(i,j-2,k,ncomp) + a(i,j-3,k,ncomp));

} else if (j == domhi[1]+1 && hi_bc[1] != Interior) {
} else if (j == domhi[1]+1 && hi_bc[1] == Symmetry) {
// use a stencil for when the interface is on the
// right physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j-1,k,ncomp) - 23.0_rt*a(i,j-2,k,ncomp) +
13.0_rt*a(i,j-3,k,ncomp) - 3.0_rt*a(i,j-4,k,ncomp));

} else if (j == domlo[1]-1 && lo_bc[1] == Outflow) {
// extrapolate to the domlo[1]-1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(i,domlo[1],k,ncomp) - 6.0_rt*a(i,domlo[1]+1,k,ncomp) + 4.0_rt*a(i,domlo[1]+2,k,ncomp) - a(i,domlo[1]+3,k,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i,j+1,k,ncomp) +
13.0_rt*a(i,j+2,k,ncomp) - 3.0_rt*a(i,j+3,k,ncomp));

} else if (j == domhi[1]+2 && hi_bc[1] == Outflow) {
// extrapolate to the domhi[1]+1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(i,domhi[1],k,ncomp) - 6.0_rt*a(i,domhi[1]-1,k,ncomp) + 4.0_rt*a(i,domhi[1]-2,k,ncomp) - a(i,domhi[1]-3,k,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i,j-2,k,ncomp) +
13.0_rt*a(i,j-3,k,ncomp) - 3.0_rt*a(i,j-4,k,ncomp));

} else {
// regular stencil
a_int(i,j,k) = (7.0_rt/12.0_rt)*(a(i,j-1,k,ncomp) + a(i,j,k,ncomp)) -
Expand All @@ -164,52 +120,30 @@ Castro::fourth_interfaces(const Box& bx,

// interpolate to the edges

if (k == domlo[2]+1 && lo_bc[2] != Interior) {
if (k == domlo[2]+1 && lo_bc[2] == Symmetry) {
// use a stencil for the interface that is one zone
// from the left physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k-1,ncomp) + 13.0_rt*a(i,j,k,ncomp) -
5.0_rt*a(i,j,k+1,ncomp) + a(i,j,k+2,ncomp));

} else if (k == domlo[2] && lo_bc[2] != Interior) {
} else if (k == domlo[2] && lo_bc[2] == Symmetry) {
// use a stencil for when the interface is on the
// left physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i,j,k+1,ncomp) +
13.0_rt*a(i,j,k+2,ncomp) - 3.0_rt*a(i,j,k+3,ncomp));

} else if (k == domhi[2] && hi_bc[2] != Interior) {
} else if (k == domhi[2] && hi_bc[2] == Symmetry) {
// use a stencil for the interface that is one zone
// from the right physical boundary, MC Eq. 22
a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i,j,k-1,ncomp) -
5.0_rt*a(i,j,k-2,ncomp) + a(i,j,k-3,ncomp));

} else if (k == domhi[2]+1 && hi_bc[2] != Interior) {
} else if (k == domhi[2]+1 && hi_bc[2] == Symmetry) {
// use a stencil for when the interface is on the
// right physical boundary MC Eq. 21
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k-1,ncomp) - 23.0_rt*a(i,j,k-2,ncomp) +
13.0_rt*a(i,j,k-3,ncomp) - 3.0_rt*a(i,j,k-4,ncomp));

} else if (k == domlo[2]-1 && lo_bc[2] == Outflow) {
// extrapolate to the domlo[2]-1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(i,j,domlo[2],ncomp) - 6.0_rt*a(i,j,domlo[2]+1,ncomp) + 4.0_rt*a(i,j,domlo[2]+2,ncomp) - a(i,j,domlo[2]+3,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i,j,k+1,ncomp) +
13.0_rt*a(i,j,k+2,ncomp) - 3.0_rt*a(i,j,k+3,ncomp));

} else if (k == domhi[2]+2 && hi_bc[2] == Outflow) {
// extrapolate to the domhi[2]+1 cell using a
// conservative cubic polynomial averaged over
// the cell
Real ac = 4.0_rt*a(i,j,domhi[2],ncomp) - 6.0_rt*a(i,j,domhi[2]-1,ncomp) + 4.0_rt*a(i,j,domhi[2]-2,ncomp) - a(i,j,domhi[2]-3,ncomp);

// now use the 1-sided stencil from above with
// this extrapolated value
a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*ac - 23.0_rt*a(i,j,k-2,ncomp) +
13.0_rt*a(i,j,k-3,ncomp) - 3.0_rt*a(i,j,k-4,ncomp));

} else {
// regular stencil
a_int(i,j,k) = (7.0_rt/12.0_rt)*(a(i,j,k-1,ncomp) + a(i,j,k,ncomp)) -
Expand Down Expand Up @@ -367,7 +301,7 @@ Castro::states(const Box& bx,
// reset the left state at domlo[0] if needed -- it is outside the domain

if (lo_bc[0] == Outflow) {
//al(domlo[0],j,k,:) = ar(domlo[0],j,k,:)
//al(domlo[0],j,k,ncomp) = ar(domlo[0],j,k,ncomp);

} else if (lo_bc[0] == Symmetry) {
if (ncomp == QU) {
Expand All @@ -391,7 +325,7 @@ Castro::states(const Box& bx,
// reset the right state at domhi[0]+1 if needed -- it is outside the domain

if (hi_bc[0] == Outflow) {
//ar(domhi[0]+1,j,k,:) = al(domhi[0]+1,j,k,:)
//ar(domhi[0]+1,j,k,ncomp) = al(domhi[0]+1,j,k,ncomp);

} else if (hi_bc[0] == Symmetry) {
if (ncomp == QU) {
Expand Down Expand Up @@ -532,7 +466,7 @@ Castro::states(const Box& bx,
// reset the left state at domlo[1] if needed -- it is outside the domain

if (lo_bc[1] == Outflow) {
//al(i,domlo[1],k,:) = ar(i,domlo[1],k,:)
//al(i,domlo[1],k,ncomp) = ar(i,domlo[1],k,ncomp);

} else if (lo_bc[1] == Symmetry) {
if (ncomp == QV) {
Expand All @@ -556,7 +490,7 @@ Castro::states(const Box& bx,
// reset the right state at domhi[1]+1 if needed -- it is outside the domain

if (hi_bc[1] == Outflow) {
//ar(i,domhi[1]+1,k,:) = al(i,domhi[1]+1,k,:)
//ar(i,domhi[1]+1,k,ncomp) = al(i,domhi[1]+1,k,ncomp);

} else if (hi_bc[1] == Symmetry) {
if (ncomp == QV) {
Expand Down Expand Up @@ -695,7 +629,7 @@ Castro::states(const Box& bx,
// reset the left state at domlo[2] if needed -- it is outside the domain

if (lo_bc[2] == Outflow) {
//al(i,j,domlo[2],:) = ar(i,j,domlo[2],:)
//al(i,j,domlo[2],ncomp) = ar(i,j,domlo[2],ncomp);

} else if (lo_bc[2] == Symmetry) {
if (ncomp == QW) {
Expand All @@ -718,7 +652,7 @@ Castro::states(const Box& bx,
// reset the right state at domhi[2]+1 if needed -- it is outside the domain

if (hi_bc[2] == Outflow) {
//ar(i,j,domhi[2]+1,:) = al(i,j,domhi[2]+1,:)
//ar(i,j,domhi[2]+1,ncomp) = al(i,j,domhi[2]+1,ncomp);

} else if (hi_bc[2] == Symmetry) {
if (ncomp == QW) {
Expand Down

0 comments on commit be6201a

Please sign in to comment.