Overset symmetry bug

[g-rydquist]

Closes #652. To copy/paste from that issue:

The overset implementation (and presumably the ALE implementation as far as I can tell?) does not properly account for the motion of the background mesh when calculating symmetry boundary conditions. As a result, the normal velocity on these boundaries of the moving mesh will set the normal velocity equal to zero, instead of the normal component of the mesh velocity. The proper way to handle this seems to be to pass the mesh velocity into the symbc() function, in a similar way to how the it is passed into, e.g., CompFlow::rhs().

The implementation is as described above: the mesh velocities are read directly into the symbc() method. The changes pushed at the time of writing do indeed affect ALECG implementation, and as a result these tests fail.

---

This change is 

[adityakpandare]

Thanks, @g-rydquist! A couple of minor comments. Also, not sure how we can get the ALECG tests to pass. I would like to keep those tests, because we might need them in the future. The ALECG tests assume that the symmetry BC is used on symmetric faces (i.e. like on a 3D mesh's y- and z-faces that is used to solve a 1D shocktube). So the zero normal-velocity (previous implementation) is enforced on those faces. It makes sense there, but not in our overset implementation. Makes me think- should we call this BC what it is- a "slip-wall" BC, and create a new function for it?

Reviewed 6 of 6 files at r1, all commit messages.
Reviewable status: all files reviewed, 3 unresolved discussions (waiting on @g-rydquist)

---

src/Inciter/NodeDiagnostics.cpp line 98 at r1 (raw file):

    // Evaluate analytic solution (if exist, if not, IC)
    auto an = u;
    auto mv = d.meshvel();

There's a bit of confusion with regards to mesh-vel that I was unaware of. Apparently there are two functions in Discretization that are called meshvel() and MeshVel(). I believe both should return the same thing for overset, but for consistency with the rest of the code, let's call the latter for Overset. Could you replace this call?

---

src/Inciter/OversetFE.cpp line 931 at r1 (raw file):

        // Apply symmetry BCs
        g_cgpde[d->MeshId()].symbc( m_u, d->meshvel(), coord, m_bnorm, m_symbcnodes );

Same here, let's use MeshVel().

---

src/Inciter/OversetFE.cpp line 1133 at r1 (raw file):

          m_triinpoel, d->Gid(), d->Bid(), d->Lid(), m_dfn, m_psup, m_esup,
          m_symbctri, d->Vol(), m_edgenode, m_edgeid,
          m_boxnodes, m_chBndGrad, m_u, d->meshvel(), m_tp, d->Boxvol(),

Same here, MeshVel.

[g-rydquist]

Yes, this definitely makes sense to me. That will be next up on my list. I replaced the meshvel()'s with the latest commits and also fixed up that bug where the coordinates were advanced from the previous RK stage instead of the coordinate location at the start of the time step.

Reviewable status: 4 of 6 files reviewed, 3 unresolved discussions (waiting on @adityakpandare)

---

src/Inciter/NodeDiagnostics.cpp line 98 at r1 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

There's a bit of confusion with regards to mesh-vel that I was unaware of. Apparently there are two functions in Discretization that are called meshvel() and MeshVel(). I believe both should return the same thing for overset, but for consistency with the rest of the code, let's call the latter for Overset. Could you replace this call?

Done.

---

src/Inciter/OversetFE.cpp line 931 at r1 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Same here, let's use MeshVel().

Done.

---

src/Inciter/OversetFE.cpp line 1133 at r1 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Same here, MeshVel.

Done.

[adityakpandare]

Thanks, @g-rydquist! I have a couple of questions about the coord updates.

Reviewed 2 of 2 files at r2.
Reviewable status: all files reviewed (commit messages unreviewed), 3 unresolved discussions (waiting on @g-rydquist)

---

src/Inciter/OversetFE.cpp line 1286 at r2 (raw file):

        // rectilinear motion
        // ---------------------------------------------------------------------

Lets move the above two comment lines below the rotatePoint() call, since that's still rotational motion

---

src/Inciter/OversetFE.cpp line 1299 at r2 (raw file):

          dsT = m_centMassVel[i]*dtp + 0.5*a_mesh[i]*dtp*dtp;
          // mesh displacement from rotation
          dsR = rCM[i] + m_centMass[i] - d->Coord()[i][p];

Not sure about this one- should rCM now be defined based on Coordn (it is currently defined based on coordK, i.e. the stage-value of the coord)?

---

src/Inciter/OversetFE.cpp line 1301 at r2 (raw file):

          dsR = rCM[i] + m_centMass[i] - d->Coord()[i][p];
          // add both contributions
          d->Coord()[i][p] = d->Coordn()[i][p] + dsT + dsR;

Great find! Now that we're using Coordn (rather than coordK) to update coord, should we also be doing the same for centMass, centMassVel, and angVelMesh so that they are consistent?

[g-rydquist]

Reviewable status: 5 of 7 files reviewed, 3 unresolved discussions (waiting on @adityakpandare and @g-rydquist)

---

src/Inciter/OversetFE.cpp line 1286 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Lets move the above two comment lines below the rotatePoint() call, since that's still rotational motion

Done.

---

src/Inciter/OversetFE.cpp line 1299 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Not sure about this one- should rCM now be defined based on Coordn (it is currently defined based on coordK, i.e. the stage-value of the coord)?

Yes I think you are absolutely correct, nice catch.

---

src/Inciter/OversetFE.cpp line 1301 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Great find! Now that we're using Coordn (rather than coordK) to update coord, should we also be doing the same for centMass, centMassVel, and angVelMesh so that they are consistent?

Yes, good catch.

[g-rydquist]

So this latest set of commits add a slipwallbc and revert the symbc to how it was previously. There's a couple things to note:

1. bndPressureint() is now only performed over slip walls, so symmetry does not contribute to the force driving the overset mesh, which I think makes sense.
2. I noticed while going through this that it seems like the fluxes calculated along the boundary with bndint() were also not accounting for the mesh movement. I did not change this at all, but I'll point out where it is in a comment below.

Reviewable status: 0 of 12 files reviewed, 5 unresolved discussions (waiting on @adityakpandare and @g-rydquist)

---

src/Inciter/ALECG.cpp line 198 at r3 (raw file):

  for (std::size_t e=0; e<m_triinpoel.size()/3; ++e)
    if (m_symbcnodes.find(m_triinpoel[e*3+0]) != end(m_symbcnodes) ||
        m_slipwallbcnodes.find(m_triinpoel[e*3+0]) != end(m_slipwallbcnodes))

This (and the statement below similar to it) I would imagine are temporary. I just wanted to get both slip and symmetry nodes into the boundary flux calculation, and didn't want to make too large of changes. But my thinking is that we should need both the slip and symmetry faces to be read into bndint() as these faces will need to be treated differently.

---

src/PDE/CompFlow/CGCompFlow.hpp line 1375 at r3 (raw file):

        p = m_mat_blk[0].compute< EOS::pressure >( rA, ruA/rA, rvA/rA, rwA/rA,
          reA );
        vn = sym ? 0.0 : (nx*(ruA/rA-w1A) + ny*(rvA/rA-w2A) + nz*(rwA/rA-w3A));

Here (and other statements like it below) are what I was referencing above. It seems as though, for slip walls, we would want to set this velocity to the normal component of the mesh velocity?