Overset symmetry bug

src/Inciter/ALECG.cpp

@@ -869,7 +869,7 @@ ALECG::BC()
       if (bc[c].first) m_u(b,c) = bc[c].second;
 
   // Apply symmetry BCs
-  g_cgpde[d->MeshId()].symbc( m_u, coord, m_bnorm, m_symbcnodes );
+  g_cgpde[d->MeshId()].symbc( m_u, d->meshvel(), coord, m_bnorm, m_symbcnodes );
 
   // Apply farfield BCs
   g_cgpde[d->MeshId()].farfieldbc( m_u, coord, m_bnorm, m_farfieldbcnodes );

src/Inciter/NodeDiagnostics.cpp

@@ -95,6 +95,7 @@ NodeDiagnostics::compute(
 
     // Evaluate analytic solution (if exist, if not, IC)
     auto an = u;
+    auto mv = d.MeshVel();
     for (std::size_t i=0; i<an.nunk(); ++i) {
       // Query analytic solution for all components of all PDEs integrated
       std::vector< tk::real > a;
@@ -104,7 +105,7 @@ NodeDiagnostics::compute(
       for (std::size_t c=0; c<an.nprop(); ++c) an(i,c) = a[c];
     }
     // Apply symmetry BCs on analytic solution (if exist, if not, IC)
-    g_cgpde[d.MeshId()].symbc( an, coord, bnorm, symbcnodes );
+    g_cgpde[d.MeshId()].symbc( an, mv, coord, bnorm, symbcnodes );
     // Apply farfield BCs on analytic solution (if exist, if not, IC)
     g_cgpde[d.MeshId()].farfieldbc( an, coord, bnorm, farfieldbcnodes );
 

src/Inciter/OversetFE.cpp

@@ -100,7 +100,10 @@ OversetFE::OversetFE( const CProxy_Discretization& disc,
   m_surfTorque({{0, 0, 0}}),
   m_centMass({{0, 0, 0}}),
   m_centMassVel({{0, 0, 0}}),
-  m_angVelMesh(0)
+  m_angVelMesh(0),
+  m_centMassn({{0, 0, 0}}),
+  m_centMassVeln({{0, 0, 0}}),
+  m_angVelMeshn(0)
 // *****************************************************************************
 //  Constructor
 //! \param[in] disc Discretization proxy
@@ -928,7 +931,7 @@ OversetFE::BC()
             if (bc[c].first) m_u(b,c) = bc[c].second;
 
         // Apply symmetry BCs
-        g_cgpde[d->MeshId()].symbc( m_u, coord, m_bnorm, m_symbcnodes );
+        g_cgpde[d->MeshId()].symbc( m_u, d->MeshVel(), coord, m_bnorm, m_symbcnodes );
 
         // Apply farfield BCs
         if (bci.get< tag::farfield >().empty() || (d->MeshId() == 0)) {
@@ -943,6 +946,17 @@ OversetFE::BC()
   }
 }
 
+void
+OversetFE::UpdateCenterOfMass()
+// *****************************************************************************
+// Update quantities associated with the center of mass at a new time step
+// *****************************************************************************
+{
+  m_centMassn = m_centMass;
+  m_centMassVeln = m_centMassVel;
+  m_angVelMeshn = m_angVelMesh;
+}
+
 void
 OversetFE::next()
 // *****************************************************************************
@@ -1196,6 +1210,8 @@ OversetFE::solve()
   // Update state at time n
   if (m_stage == 0) {
     m_un = m_u;
+    d->UpdateCoordn();
+    UpdateCenterOfMass();
   }
 
   // Explicit time-stepping using RK3
@@ -1257,9 +1273,9 @@ OversetFE::solve()
         // rotation (this is currently only configured for planar motion)
         // ---------------------------------------------------------------------
         std::array< tk::real, 3 > rCM{{
-          d->Coord()[0][p] - m_centMass[0],
-          d->Coord()[1][p] - m_centMass[1],
-          d->Coord()[2][p] - m_centMass[2] }};
+          d->Coord()[0][p] - m_centMassn[0],
+          d->Coord()[1][p] - m_centMassn[1],
+          d->Coord()[2][p] - m_centMassn[2] }};
 
         // obtain tangential velocity
         tk::real r_mag(0.0);
@@ -1281,39 +1297,39 @@ OversetFE::solve()
         // angle of rotation
         auto dtheta = m_angVelMesh*dtp + 0.5*alpha_mesh*dtp*dtp;
 
+        // add contribution of rotation to mesh displacement
+        std::array< tk::real, 3 > angles{{ 0, 0, 0 }};
+        angles[sym_dir] = dtheta * 180.0/pi;
+        tk::rotatePoint(angles, rCM);
+
         // rectilinear motion
         // ---------------------------------------------------------------------
-        std::array< tk::real, 3 > ds{{ 0, 0, 0 }};
+
+        tk::real dsT, dsR;
+
         for (std::size_t i=0; i<3; ++i) {
-          // mesh displacement
-          ds[i] = m_centMassVel[i]*dtp + 0.5*a_mesh[i]*dtp*dtp;
-          // mesh velocity
+          // mesh displacement from translation
+          dsT = m_centMassVel[i]*dtp + 0.5*a_mesh[i]*dtp*dtp;
+          // mesh displacement from rotation
+          dsR = rCM[i] + m_centMass[i] - d->Coordn()[i][p];
+          // add both contributions
+          d->Coord()[i][p] = d->Coordn()[i][p] + dsT + dsR;
+          // mesh velocity change from translation
           u_mesh(p,i) += a_mesh[i]*dtp;
         }
 
         // add contribution of rotation to mesh velocity
         u_mesh(p,i1) += a1*dtp;
         u_mesh(p,i2) += a2*dtp;
-
-        // add contribution of rotation to mesh displacement
-        std::array< tk::real, 3 > angles{{ 0, 0, 0 }};
-        angles[sym_dir] = dtheta * 180.0/pi;
-        tk::rotatePoint(angles, rCM);
-
-        // combine both contributions
-        for (std::size_t i=0; i<3; ++i) {
-          d->Coord()[i][p] = rCM[i] + m_centMass[i];
-          d->Coord()[i][p] += ds[i];
-        }
       }
 
       // update angular velocity
-      m_angVelMesh += alpha_mesh*dtp;
+      m_angVelMesh = m_angVelMeshn + alpha_mesh*dtp;
 
       // move center of mass
       for (std::size_t i=0; i<3; ++i) {
-        m_centMass[i] += m_centMassVel[i]*dtp + 0.5*a_mesh[i]*dtp*dtp;
-        m_centMassVel[i] += a_mesh[i]*dtp;  // no rotational component
+        m_centMass[i] = m_centMassn[i] + m_centMassVel[i]*dtp + 0.5*a_mesh[i]*dtp*dtp;
+        m_centMassVel[i] = m_centMassVeln[i] + a_mesh[i]*dtp;  // no rotational component
       }
     }
 

src/Inciter/OversetFE.hpp

@@ -242,6 +242,9 @@ class OversetFE : public CBase_OversetFE {
       p | m_centMass;
       p | m_centMassVel;
       p | m_angVelMesh;
+      p | m_centMassn;
+      p | m_centMassVeln;
+      p | m_angVelMeshn;
     }
     //! \brief Pack/Unpack serialize operator|
     //! \param[in,out] p Charm++'s PUP::er serializer object reference
@@ -381,6 +384,12 @@ class OversetFE : public CBase_OversetFE {
     std::array< tk::real, 3 > m_centMassVel;
     //! Angular velocity of the rigid body
     tk::real m_angVelMesh;
+    //! Center of mass of rigid body at time n
+    std::array< tk::real, 3 > m_centMassn;
+    //! Velocity of the center of mass of rigid body at time n
+    std::array< tk::real, 3 > m_centMassVeln;
+    //! Angular velocity of the rigid body at time n
+    tk::real m_angVelMeshn;
 
     //! Access bound Discretization class pointer
     Discretization* Disc() const {
@@ -453,6 +462,9 @@ class OversetFE : public CBase_OversetFE {
 
     //! Apply boundary conditions
     void BC();
+
+    //! Update quantities associated with the center of mass at a new time step
+    void UpdateCenterOfMass();
 };
 
 } // inciter::

src/PDE/CGPDE.hpp

@@ -218,12 +218,13 @@ class CGPDE {
     //! Public interface to set symmetry boundary conditions at nodes
     void
     symbc( tk::Fields& U,
+           const tk::Fields& W,
            const std::array< std::vector< real >, 3 >& coord,
            const std::unordered_map< int,
                    std::unordered_map< std::size_t,
                      std::array< real, 4 > > >& bnorm,
            const std::unordered_set< std::size_t >& nodes ) const
-    { self->symbc( U, coord, bnorm, nodes ); }
+    { self->symbc( U, W, coord, bnorm, nodes ); }
 
     //! Public interface to set farfield boundary conditions at nodes
     void
@@ -385,6 +386,7 @@ class CGPDE {
              bool ) const = 0;
       virtual void symbc(
         tk::Fields& U,
+        const tk::Fields& W,
         const std::array< std::vector< real >, 3 >&,
         const std::unordered_map< int,
                 std::unordered_map< std::size_t,
@@ -520,12 +522,13 @@ class CGPDE {
                                       increment ); }
       void symbc(
         tk::Fields& U,
+        const tk::Fields& W,
         const std::array< std::vector< real >, 3 >& coord,
         const std::unordered_map< int,
                 std::unordered_map< std::size_t,
                   std::array< real, 4 > > >& bnorm,
         const std::unordered_set< std::size_t >& nodes ) const override
-      { data.symbc( U, coord, bnorm, nodes ); }
+      { data.symbc( U, W, coord, bnorm, nodes ); }
       void farfieldbc(
         tk::Fields& U,
         const std::array< std::vector< real >, 3 >& coord,

src/PDE/CompFlow/CGCompFlow.hpp

@@ -750,11 +750,13 @@ class CompFlow {
 
     //! Set symmetry boundary conditions at nodes
     //! \param[in] U Solution vector at recent time step
+    //! \param[in] W Mesh velocity
     //! \param[in] bnorm Face normals in boundary points, key local node id,
     //!   first 3 reals of value: unit normal, outer key: side set id
     //! \param[in] nodes Unique set of node ids at which to set symmetry BCs
     void
     symbc( tk::Fields& U,
+           const tk::Fields& W,
            const std::array< std::vector< real >, 3 >&,
            const std::unordered_map< int,
              std::unordered_map< std::size_t, std::array< real, 4 > > >& bnorm,
@@ -778,12 +780,12 @@ class CompFlow {
               if (i != end(j->second)) {
                 std::array< real, 3 >
                   n{ i->second[0], i->second[1], i->second[2] },
-                  v{ U(p,1), U(p,2), U(p,3) };
-                auto v_dot_n = tk::dot( v, n );
-                // symbc: remove normal component of velocity
-                U(p,1) -= v_dot_n * n[0];
-                U(p,2) -= v_dot_n * n[1];
-                U(p,3) -= v_dot_n * n[2];
+                  rel_v{ U(p,1) - W(p,0), U(p,2) - W(p,1), U(p,3) - W(p,2) };
+                auto rel_v_dot_n = tk::dot( rel_v, n );
+                // symbc: remove normal component of relative velocity
+                U(p,1) -= rel_v_dot_n * n[0];
+                U(p,2) -= rel_v_dot_n * n[1];
+                U(p,3) -= rel_v_dot_n * n[2];
               }
             }
           }

src/PDE/Transport/CGTransport.hpp

@@ -398,6 +398,7 @@ class Transport {
     void
     symbc(
       tk::Fields&,
+      const tk::Fields&,
       const std::array< std::vector< real >, 3 >&,
       const std::unordered_map< int,
               std::unordered_map< std::size_t,