This is not fully clear to me, and if we clobber this solver with the semi-Lagrangian approach it may not be worth looking into too deeply.
The fluxes in AdvDiffSystem are stored as a flat vector representing the 2D grid in column-major order. In AdvDiffSystem::forwardEulerAdvection we loop over all indexes in this flat vector, to compute the 2nd order flux. For each index, we check if it is on a boundary (special treatment), then find the relevant direction for the upwind scheme (u_local or v_local check) and apply the 3-point upwind stencil.
The problem lies in the minmod_X_YPos series of functions. When they are called, we've checked that this edge is not a boundary. However because we use a 3-point stencil, an edge that is 1-away from the boundary can still reach outside of the grid when computing it's flux. The point that's outside the boundary should be a ghost node. However the call to minmod_N_vNeg for example does not realize this, and reaches for these points for stencil:
double phi_P = phi_[pointID];
double phi_N = phi_[pointID + 1];
double phi_NN = phi_[pointID + 2];
but here phi_NN is the bottom of the next column instead of a ghost node. I think this occurs for (almost?) all edge by reading the adjacent point in memory but that is not the adjacent point on the grid.
This only affects the edges, but could maybe explain some of the weird behaviors we have seen with the contrail mask? This essentially teleports mass from one part of the grid to an edge, and if it gets bad enough maybe this would slowly expand the mask.
This is not fully clear to me, and if we clobber this solver with the semi-Lagrangian approach it may not be worth looking into too deeply.
The fluxes in
AdvDiffSystemare stored as a flat vector representing the 2D grid in column-major order. InAdvDiffSystem::forwardEulerAdvectionwe loop over all indexes in this flat vector, to compute the 2nd order flux. For each index, we check if it is on a boundary (special treatment), then find the relevant direction for the upwind scheme (u_localorv_localcheck) and apply the 3-point upwind stencil.The problem lies in the
minmod_X_YPosseries of functions. When they are called, we've checked that this edge is not a boundary. However because we use a 3-point stencil, an edge that is 1-away from the boundary can still reach outside of the grid when computing it's flux. The point that's outside the boundary should be a ghost node. However the call tominmod_N_vNegfor example does not realize this, and reaches for these points for stencil:but here
phi_NNis the bottom of the next column instead of a ghost node. I think this occurs for (almost?) all edge by reading the adjacent point in memory but that is not the adjacent point on the grid.This only affects the edges, but could maybe explain some of the weird behaviors we have seen with the contrail mask? This essentially teleports mass from one part of the grid to an edge, and if it gets bad enough maybe this would slowly expand the mask.