Skip to content

Commit

Permalink
Merge pull request #166 from Pressio/fix_bc_stencil_weno
Browse files Browse the repository at this point in the history 
fix ghost values for weno3,5, pointed out by Chris
  • Loading branch information
@fnrizzi
fnrizzi committed Dec 20, 2023
2 parents 6a6f4b1 + cb5d4b8 commit fb15928
Show file tree
Hide file tree
Showing 60 changed files with 15,240 additions and 13,261 deletions.
256 changes: 143 additions & 113 deletions include/pressiodemoapps/impl/euler_2d_ghost_filler_neumann.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,15 +79,27 @@ class Ghost2dNeumannFiller
template<class index_t>
void operator()(index_t smPt, int gRow)
{
if (m_numDofPerCell == 4){
fourDofImpl(smPt, gRow);
if (m_stencilSize == 3){
stencilThreeImpl(smPt, gRow);
}

else if (m_stencilSize == 5){
stencilFiveImpl(smPt, gRow);
}

else if (m_stencilSize == 7){
stencilSevenImpl(smPt, gRow);
}

else{
throw std::runtime_error("euler riemann2d ghost filler: invalid stencil size");
}
}

private:

template<class index_t>
void fourDofImpl(index_t smPt, int gRow)
void stencilThreeImpl(index_t smPt, int gRow)
{
const auto & graph = m_meshObj.graph();
const auto cellGID = graph(smPt, 0);
Expand Down Expand Up @@ -129,119 +141,137 @@ class Ghost2dNeumannFiller
m_ghostBack(gRow, 2) = m_state(uIndex+2);
m_ghostBack(gRow, 3) = m_state(uIndex+3);
}
}

template<class index_t>
void stencilFiveImpl(index_t smPt, int gRow)
{
const auto & graph = m_meshObj.graph();
const auto cellGID = graph(smPt, 0);
const auto uIndex = cellGID*m_numDofPerCell;

stencilThreeImpl(smPt, gRow);
const auto left0 = graph(smPt, 1);
const auto front0 = graph(smPt, 2);
const auto right0 = graph(smPt, 3);
const auto back0 = graph(smPt, 4);
const auto left1 = graph(smPt, 5);
const auto front1 = graph(smPt, 6);
const auto right1 = graph(smPt, 7);
const auto back1 = graph(smPt, 8);

if (left1 == -1){
auto ind = uIndex;
if (left0==-1){ ind = right0*m_numDofPerCell; }
else { ind = left0*m_numDofPerCell; }

m_ghostLeft(gRow, 4) = m_state(ind);
m_ghostLeft(gRow, 5) = m_state(ind+1);
m_ghostLeft(gRow, 6) = m_state(ind+2);
m_ghostLeft(gRow, 7) = m_state(ind+3);
}

if (front1 == -1){
auto ind = uIndex;
if (front0==-1){ ind = back0*m_numDofPerCell; }
else { ind = front0*m_numDofPerCell; }

m_ghostFront(gRow, 4) = m_state(ind);
m_ghostFront(gRow, 5) = m_state(ind+1);
m_ghostFront(gRow, 6) = m_state(ind+2);
m_ghostFront(gRow, 7) = m_state(ind+3);
}

if (right1 == -1){
auto ind = uIndex;
if (right0==-1){ ind = left0*m_numDofPerCell; }
else { ind = right0*m_numDofPerCell; }

m_ghostRight(gRow, 4) = m_state(ind);
m_ghostRight(gRow, 5) = m_state(ind+1);
m_ghostRight(gRow, 6) = m_state(ind+2);
m_ghostRight(gRow, 7) = m_state(ind+3);
}

if (m_stencilSize >= 5){
const auto left1 = graph(smPt, 5);
const auto front1 = graph(smPt, 6);
const auto right1 = graph(smPt, 7);
const auto back1 = graph(smPt, 8);

if (left1 == -1){
const auto ind = right0*m_numDofPerCell;
m_ghostLeft(gRow, 4) = m_state(ind);
m_ghostLeft(gRow, 5) = m_state(ind+1);
m_ghostLeft(gRow, 6) = m_state(ind+2);
m_ghostLeft(gRow, 7) = m_state(ind+3);
}

if (front1 == -1){
const auto ind = back0*m_numDofPerCell;
m_ghostFront(gRow, 4) = m_state(ind);
m_ghostFront(gRow, 5) = m_state(ind+1);
m_ghostFront(gRow, 6) = m_state(ind+2);
m_ghostFront(gRow, 7) = m_state(ind+3);
}

if (right1 == -1){
const auto ind = left0*m_numDofPerCell;
m_ghostRight(gRow, 4) = m_state(ind);
m_ghostRight(gRow, 5) = m_state(ind+1);
m_ghostRight(gRow, 6) = m_state(ind+2);
m_ghostRight(gRow, 7) = m_state(ind+3);
}

if (back1 == -1){
const auto ind = front0*m_numDofPerCell;
m_ghostBack(gRow, 4) = m_state(ind);
m_ghostBack(gRow, 5) = m_state(ind+1);
m_ghostBack(gRow, 6) = m_state(ind+2);
m_ghostBack(gRow, 7) = m_state(ind+3);
}
if (back1 == -1){
auto ind = uIndex;
if (back0==-1){ ind = front0*m_numDofPerCell; }
else { ind = back0*m_numDofPerCell; }

m_ghostBack(gRow, 4) = m_state(ind);
m_ghostBack(gRow, 5) = m_state(ind+1);
m_ghostBack(gRow, 6) = m_state(ind+2);
m_ghostBack(gRow, 7) = m_state(ind+3);
}
}

template<class index_t>
void stencilSevenImpl(index_t smPt, int gRow)
{
const auto & graph = m_meshObj.graph();
const auto cellGID = graph(smPt, 0);
const auto uIndex = cellGID*m_numDofPerCell;

stencilFiveImpl(smPt, gRow);
const auto left0 = graph(smPt, 1);
const auto front0 = graph(smPt, 2);
const auto right0 = graph(smPt, 3);
const auto back0 = graph(smPt, 4);
const auto left1 = graph(smPt, 5);
const auto front1 = graph(smPt, 6);
const auto right1 = graph(smPt, 7);
const auto back1 = graph(smPt, 8);
const auto left2 = graph(smPt, 9);
const auto front2 = graph(smPt, 10);
const auto right2 = graph(smPt, 11);
const auto back2 = graph(smPt, 12);

if (left2 == -1){
auto ind = uIndex; ;
if (left1!=-1 && left0!=-1){ ind = left1*m_numDofPerCell; }
if (left1==-1 && left0!=-1){ ind = uIndex; }
if (left1==-1 && left0==-1){ ind = right1*m_numDofPerCell; }

m_ghostLeft(gRow, 8) = m_state(ind);
m_ghostLeft(gRow, 9) = m_state(ind+1);
m_ghostLeft(gRow, 10) = m_state(ind+2);
m_ghostLeft(gRow, 11) = m_state(ind+3);
}

if (front2 == -1){
auto ind = uIndex; ;
if (front1!=-1 && front0!=-1){ ind = front1*m_numDofPerCell; }
if (front1==-1 && front0!=-1){ ind = uIndex; }
if (front1==-1 && front0==-1){ ind = back1*m_numDofPerCell; }

m_ghostFront(gRow, 8) = m_state(ind);
m_ghostFront(gRow, 9) = m_state(ind+1);
m_ghostFront(gRow, 10) = m_state(ind+2);
m_ghostFront(gRow, 11) = m_state(ind+3);
}

if (right2 == -1){
auto ind = uIndex; ;
if (right1!=-1 && right0!=-1){ ind = right1*m_numDofPerCell; }
if (right1==-1 && right0!=-1){ ind = uIndex; }
if (right1==-1 && right0==-1){ ind = left1*m_numDofPerCell; }

m_ghostRight(gRow, 8) = m_state(ind);
m_ghostRight(gRow, 9) = m_state(ind+1);
m_ghostRight(gRow, 10) = m_state(ind+2);
m_ghostRight(gRow, 11) = m_state(ind+3);
}

if (back2 == -1){
auto ind = uIndex; ;
if (back1!=-1 && back0!=-1){ ind = back1*m_numDofPerCell; }
if (back1==-1 && back0!=-1){ ind = uIndex; }
if (back1==-1 && back0==-1){ ind = front1*m_numDofPerCell; }

if (m_stencilSize == 7){
const auto left1 = graph(smPt, 5);
const auto front1 = graph(smPt, 6);
const auto right1 = graph(smPt, 7);
const auto back1 = graph(smPt, 8);
const auto left2 = graph(smPt, 9);
const auto front2 = graph(smPt, 10);
const auto right2 = graph(smPt, 11);
const auto back2 = graph(smPt, 12);

if (left1 == -1){
const auto ind = right0*m_numDofPerCell;
m_ghostLeft(gRow, 4) = m_state(ind);
m_ghostLeft(gRow, 5) = m_state(ind+1);
m_ghostLeft(gRow, 6) = m_state(ind+2);
m_ghostLeft(gRow, 7) = m_state(ind+3);
}

if (front1 == -1){
const auto ind = back0*m_numDofPerCell;
m_ghostFront(gRow, 4) = m_state(ind);
m_ghostFront(gRow, 5) = m_state(ind+1);
m_ghostFront(gRow, 6) = m_state(ind+2);
m_ghostFront(gRow, 7) = m_state(ind+3);
}

if (right1 == -1){
const auto ind = left0*m_numDofPerCell;
m_ghostRight(gRow, 4) = m_state(ind);
m_ghostRight(gRow, 5) = m_state(ind+1);
m_ghostRight(gRow, 6) = m_state(ind+2);
m_ghostRight(gRow, 7) = m_state(ind+3);
}

if (back1 == -1){
const auto ind = front0*m_numDofPerCell;
m_ghostBack(gRow, 4) = m_state(ind);
m_ghostBack(gRow, 5) = m_state(ind+1);
m_ghostBack(gRow, 6) = m_state(ind+2);
m_ghostBack(gRow, 7) = m_state(ind+3);
}

if (left2 == -1){
const auto ind = right1*m_numDofPerCell;
m_ghostLeft(gRow, 8) = m_state(ind);
m_ghostLeft(gRow, 9) = m_state(ind+1);
m_ghostLeft(gRow, 10) = m_state(ind+2);
m_ghostLeft(gRow, 11) = m_state(ind+3);
}

if (front2 == -1){
const auto ind = back1*m_numDofPerCell;
m_ghostFront(gRow, 8) = m_state(ind);
m_ghostFront(gRow, 9) = m_state(ind+1);
m_ghostFront(gRow, 10) = m_state(ind+2);
m_ghostFront(gRow, 11) = m_state(ind+3);
}

if (right2 == -1){
const auto ind = left1*m_numDofPerCell;
m_ghostRight(gRow, 8) = m_state(ind);
m_ghostRight(gRow, 9) = m_state(ind+1);
m_ghostRight(gRow, 10) = m_state(ind+2);
m_ghostRight(gRow, 11) = m_state(ind+3);
}

if (back2 == -1){
const auto ind = front1*m_numDofPerCell;
m_ghostBack(gRow, 8) = m_state(ind);
m_ghostBack(gRow, 9) = m_state(ind+1);
m_ghostBack(gRow, 10) = m_state(ind+2);
m_ghostBack(gRow, 11) = m_state(ind+3);
}
m_ghostBack(gRow, 8) = m_state(ind);
m_ghostBack(gRow, 9) = m_state(ind+1);
m_ghostBack(gRow, 10) = m_state(ind+2);
m_ghostBack(gRow, 11) = m_state(ind+3);
}
}

Expand Down
15 changes: 14 additions & 1 deletion include/pressiodemoapps/impl/euler_2d_prob_class.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,13 @@ class EigenApp
}
}

#ifdef PRESSIODEMOAPPS_ENABLE_TESTS
const auto & viewGhostLeft() const{ return m_ghostLeft; }
const auto & viewGhostFront() const{ return m_ghostFront; }
const auto & viewGhostRight() const{ return m_ghostRight; }
const auto & viewGhostBack() const{ return m_ghostBack; }
#endif

public:
template <class T>
void setBCPointer(::pressiodemoapps::impl::GhostRelativeLocation rloc, T* ptr) {
Expand Down Expand Up @@ -474,7 +481,7 @@ class EigenApp
#pragma omp for schedule(static)
#endif
for (decltype(rowsBd.size()) it=0; it<rowsBd.size(); ++it){
ghF(rowsBd[it], it);
ghF(rowsBd[it], it);
}
}

Expand Down Expand Up @@ -1255,6 +1262,12 @@ class EigenApp
::pressiodemoapps::resize(m_ghostFront, s1, numGhostValues);
::pressiodemoapps::resize(m_ghostRight, s1, numGhostValues);
::pressiodemoapps::resize(m_ghostBack, s1, numGhostValues);

constexpr auto val = std::numeric_limits<scalar_type>::min();
m_ghostLeft.setConstant(val);
m_ghostFront.setConstant(val);
m_ghostRight.setConstant(val);
m_ghostBack.setConstant(val);
}

protected:
Expand Down
Loading

0 comments on commit fb15928

Please sign in to comment.