diff --git a/src/Control/Inciter/InputDeck/InputDeck.hpp b/src/Control/Inciter/InputDeck/InputDeck.hpp index 4551edca55..dea62d60de 100644 --- a/src/Control/Inciter/InputDeck/InputDeck.hpp +++ b/src/Control/Inciter/InputDeck/InputDeck.hpp @@ -50,7 +50,8 @@ using bclist = tk::TaggedTuple< brigand::list< tag::outlet, std::vector< std::size_t >, tag::farfield, std::vector< std::size_t >, tag::extrapolate, std::vector< std::size_t >, - tag::noslipwall, std::vector< std::size_t > + tag::noslipwall, std::vector< std::size_t >, + tag::slipwall, std::vector< std::size_t > > >; // Transport @@ -146,6 +147,7 @@ using bcList = tk::TaggedTuple< brigand::list< tag::farfield, std::vector< std::size_t >, tag::extrapolate, std::vector< std::size_t >, tag::noslipwall, std::vector< std::size_t >, + tag::slipwall, std::vector< std::size_t >, tag::velocity, std::vector< tk::real >, tag::pressure, tk::real, tag::density, tk::real, @@ -1734,6 +1736,11 @@ class InputDeck : public tk::TaggedTuple< ConfigMembers > { R"(This keyword is used to list (multiple) no-slip wall BC sidesets.)", "vector of uint(s)"}); + keywords.insert({"slipwall", + "List sidesets with slip wall boundary conditions", + R"(This keyword is used to list (multiple) slip wall BC sidesets.)", + "vector of uint(s)"}); + keywords.insert({"timedep", "Start configuration block describing time dependent boundary conditions", R"(This keyword is used to introduce a bc_timedep block, used to diff --git a/src/Control/Inciter/InputDeck/LuaParser.cpp b/src/Control/Inciter/InputDeck/LuaParser.cpp index d738513428..ba0af08ad4 100644 --- a/src/Control/Inciter/InputDeck/LuaParser.cpp +++ b/src/Control/Inciter/InputDeck/LuaParser.cpp @@ -1253,6 +1253,9 @@ LuaParser::storeInputDeck( storeVecIfSpecd< uint64_t >(sol_bc[i+1], "noslipwall", bc_deck[i].get< tag::noslipwall >(), {}); + storeVecIfSpecd< uint64_t >(sol_bc[i+1], "slipwall", + bc_deck[i].get< tag::slipwall >(), {}); + // Time-dependent BC if (sol_bc[i+1]["timedep"].valid()) { const sol::table& sol_tdbc = sol_bc[i+1]["timedep"]; diff --git a/src/Control/Tags.hpp b/src/Control/Tags.hpp index f231ff0ea5..802ae4d8f5 100644 --- a/src/Control/Tags.hpp +++ b/src/Control/Tags.hpp @@ -291,6 +291,7 @@ DEFTAG(farfield); DEFTAG(extrapolate); DEFTAG(noslipwall); DEFTAG(timedep); +DEFTAG(slipwall); DEFTAG(back_pressure); DEFTAG(rigid_body_motion); diff --git a/src/Inciter/ALECG.cpp b/src/Inciter/ALECG.cpp index aa53e57aa9..cbd142ad48 100644 --- a/src/Inciter/ALECG.cpp +++ b/src/Inciter/ALECG.cpp @@ -81,7 +81,9 @@ ALECG::ALECG( const CProxy_Discretization& disc, m_bnormc(), m_symbcnodes(), m_farfieldbcnodes(), + m_slipwallbcnodes(), m_symbctri(), + m_slipwallbctri(), m_timedepbcnodes(), m_timedepbcFn(), m_stage( 0 ), @@ -173,15 +175,36 @@ ALECG::queryBnd() for (const auto& [s,nodes] : far) m_farfieldbcnodes.insert( begin(nodes), end(nodes) ); - // If farfield BC is set on a node, will not also set symmetry BC - for (auto fn : m_farfieldbcnodes) m_symbcnodes.erase(fn); + // Prepare unique set of slip wall BC nodes + auto slip = d->bcnodes< tag::slipwall >( m_bface, m_triinpoel ); + for (const auto& [s,nodes] : slip) + m_slipwallbcnodes.insert( begin(nodes), end(nodes) ); - // Prepare boundary nodes contiguously accessible from a triangle-face loop + // If farfield BC is set on a node, will not also set symmetry and slip BC + for (auto fn : m_farfieldbcnodes) { + m_symbcnodes.erase(fn); + m_slipwallbcnodes.erase(fn); + } + + // If symmetry BC is set on a node, will not also set slip BC + for (auto fn : m_symbcnodes) { + m_slipwallbcnodes.erase(fn); + } + + // Prepare boundary nodes contiguously accessible from a triangle-face loop, + // which contain both symmetry and no slip walls m_symbctri.resize( m_triinpoel.size()/3, 0 ); for (std::size_t e=0; eMeshId()].farfieldbc( m_u, coord, m_bnorm, m_farfieldbcnodes ); + // Apply slip wall BCs + g_cgpde[d->MeshId()].slipwallbc( m_u, d->meshvel(), coord, m_bnorm, + m_slipwallbcnodes ); + // Apply user defined time dependent BCs g_cgpde[d->MeshId()].timedepbc( d->T(), m_u, m_timedepbcnodes, m_timedepbcFn ); @@ -1044,7 +1071,7 @@ ALECG::rhs() conserved( m_u, Disc()->Vol() ); g_cgpde[d->MeshId()].rhs( d->T() + prev_rkcoef * d->Dt(), d->Coord(), d->Inpoel(), m_triinpoel, d->Gid(), d->Bid(), d->Lid(), m_dfn, m_psup, m_esup, - m_symbctri, d->Vol(), m_edgenode, m_edgeid, + m_symbctri, m_slipwallbctri, d->Vol(), m_edgenode, m_edgeid, m_boxnodes, m_chBndGrad, m_u, d->meshvel(), m_tp, d->Boxvol(), m_rhs ); volumetric( m_u, Disc()->Vol() ); @@ -1201,7 +1228,8 @@ ALECG::ale() conserved( m_u, Disc()->Vol() ); conserved( m_un, Disc()->Voln() ); auto diag_computed = m_diag.compute( *d, m_u, m_un, m_bnorm, - m_symbcnodes, m_farfieldbcnodes ); + m_symbcnodes, m_farfieldbcnodes, + m_slipwallbcnodes ); volumetric( m_u, Disc()->Vol() ); volumetric( m_un, Disc()->Voln() ); // Increase number of iterations and physical time diff --git a/src/Inciter/ALECG.hpp b/src/Inciter/ALECG.hpp index 470f9e58b1..a556b35e98 100644 --- a/src/Inciter/ALECG.hpp +++ b/src/Inciter/ALECG.hpp @@ -234,7 +234,9 @@ class ALECG : public CBase_ALECG { p | m_bnormc; p | m_symbcnodes; p | m_farfieldbcnodes; + p | m_slipwallbcnodes; p | m_symbctri; + p | m_slipwallbctri; p | m_timedepbcnodes; p | m_timedepbcFn; p | m_stage; @@ -334,8 +336,12 @@ class ALECG : public CBase_ALECG { std::unordered_set< std::size_t > m_symbcnodes; //! Unique set of nodes at which farfield BCs are set std::unordered_set< std::size_t > m_farfieldbcnodes; + //! Unique set of nodes at which slip wall BCs are set + std::unordered_set< std::size_t > m_slipwallbcnodes; //! Vector with 1 at symmetry BC boundary triangles std::vector< int > m_symbctri; + //! Vector with 1 at slip wall BC boundary triangles + std::vector< int > m_slipwallbctri; //! \brief Unique set of nodes at which time dependent BCs are set // for each time dependent BC std::vector< std::unordered_set< std::size_t > > m_timedepbcnodes; diff --git a/src/Inciter/NodeDiagnostics.cpp b/src/Inciter/NodeDiagnostics.cpp index 393465011a..d473b566c9 100644 --- a/src/Inciter/NodeDiagnostics.cpp +++ b/src/Inciter/NodeDiagnostics.cpp @@ -53,7 +53,8 @@ NodeDiagnostics::compute( const std::unordered_map< int, std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm, const std::unordered_set< std::size_t >& symbcnodes, - const std::unordered_set< std::size_t >& farfieldbcnodes ) const + const std::unordered_set< std::size_t >& farfieldbcnodes, + const std::unordered_set< std::size_t >& slipwallbcnodes ) const // ***************************************************************************** // Compute diagnostics, e.g., residuals, norms of errors, etc. //! \param[in] d Discretization proxy to read from @@ -64,6 +65,7 @@ NodeDiagnostics::compute( //! \param[in] symbcnodes Unique set of node ids at which to set symmetry BCs //! \param[in] farfieldbcnodes Unique set of node ids at which to set farfield //! BCs +//! \param[in] slipwallbcnodes Unique set of node ids at which to set slip BCs //! \return True if diagnostics have been computed //! \details Diagnostics are defined as some norm, e.g., L2 norm, of a quantity, //! computed in mesh nodes, A, as ||A||_2 = sqrt[ sum_i(A_i)^2 V_i ], @@ -95,6 +97,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 a; @@ -107,6 +110,8 @@ NodeDiagnostics::compute( g_cgpde[d.MeshId()].symbc( an, coord, bnorm, symbcnodes ); // Apply farfield BCs on analytic solution (if exist, if not, IC) g_cgpde[d.MeshId()].farfieldbc( an, coord, bnorm, farfieldbcnodes ); + // Apply slip wall BCs on analytic solution (if exist, if not, IC) + g_cgpde[d.MeshId()].slipwallbc( an, mv, coord, bnorm, slipwallbcnodes ); // Put in norms sweeping our mesh chunk for (std::size_t i=0; i > >& bnorm, const std::unordered_set< std::size_t >& symbcnodes, - const std::unordered_set< std::size_t >& farfieldbcnodes ) const; + const std::unordered_set< std::size_t >& farfieldbcnodes, + const std::unordered_set< std::size_t >& slipwallbcnodes ) const; /** @name Charm++ pack/unpack serializer member functions */ ///@{ diff --git a/src/Inciter/OversetFE.cpp b/src/Inciter/OversetFE.cpp index 018df9f339..b0ebfb4788 100644 --- a/src/Inciter/OversetFE.cpp +++ b/src/Inciter/OversetFE.cpp @@ -81,7 +81,9 @@ OversetFE::OversetFE( const CProxy_Discretization& disc, m_bnormc(), m_symbcnodes(), m_farfieldbcnodes(), + m_slipwallbcnodes(), m_symbctri(), + m_slipwallbctri(), m_timedepbcnodes(), m_timedepbcFn(), m_stage( 0 ), @@ -100,7 +102,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 @@ -178,15 +183,36 @@ OversetFE::getBCNodes() for (const auto& [s,nodes] : far) m_farfieldbcnodes.insert( begin(nodes), end(nodes) ); - // If farfield BC is set on a node, will not also set symmetry BC - for (auto fn : m_farfieldbcnodes) m_symbcnodes.erase(fn); + // Prepare unique set of slip wall BC nodes + auto slip = d->bcnodes< tag::slipwall >( m_bface, m_triinpoel ); + for (const auto& [s,nodes] : slip) + m_slipwallbcnodes.insert( begin(nodes), end(nodes) ); - // Prepare boundary nodes contiguously accessible from a triangle-face loop + // If farfield BC is set on a node, will not also set symmetry and slip BC + for (auto fn : m_farfieldbcnodes) { + m_symbcnodes.erase(fn); + m_slipwallbcnodes.erase(fn); + } + + // If symmetry BC is set on a node, will not also set slip BC + for (auto fn : m_symbcnodes) { + m_slipwallbcnodes.erase(fn); + } + + // Prepare boundary nodes contiguously accessible from a triangle-face loop, + // which contain both symmetry and no slip walls m_symbctri.resize( m_triinpoel.size()/3, 0 ); for (std::size_t e=0; eMeshId()].farfieldbc( m_u, coord, m_bnorm, m_farfieldbcnodes ); } + // Apply slip wall BCs + g_cgpde[d->MeshId()].slipwallbc( m_u, d->MeshVel(), coord, m_bnorm, + m_slipwallbcnodes ); + // Apply user defined time dependent BCs g_cgpde[d->MeshId()].timedepbc( d->T(), m_u, m_timedepbcnodes, m_timedepbcFn ); @@ -943,6 +973,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() // ***************************************************************************** @@ -1002,7 +1043,7 @@ OversetFE::dt() std::vector< tk::real > F(6, 0.0); if (g_inputdeck.get< tag::rigid_body_motion >().get< tag::rigid_body_movt >() && d->MeshId() > 0) { - g_cgpde[d->MeshId()].bndPressureInt( d->Coord(), m_triinpoel, m_symbctri, + g_cgpde[d->MeshId()].bndPressureInt( d->Coord(), m_triinpoel, m_slipwallbctri, m_u, m_centMass, F ); } @@ -1129,7 +1170,7 @@ OversetFE::rhs() for (std::size_t p=0; pMeshId()].rhs( d->T() + prev_rkcoef * d->Dt(), d->Coord(), d->Inpoel(), m_triinpoel, d->Gid(), d->Bid(), d->Lid(), m_dfn, m_psup, m_esup, - m_symbctri, d->Vol(), m_edgenode, m_edgeid, + m_symbctri, m_slipwallbctri, d->Vol(), m_edgenode, m_edgeid, m_boxnodes, m_chBndGrad, m_u, d->MeshVel(), m_tp, d->Boxvol(), m_rhs ); if (steady) @@ -1196,6 +1237,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 +1300,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->Coordn()[0][p] - m_centMassn[0], + d->Coordn()[1][p] - m_centMassn[1], + d->Coordn()[2][p] - m_centMassn[2] }}; // obtain tangential velocity tk::real r_mag(0.0); @@ -1281,39 +1324,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 } } @@ -1334,7 +1377,8 @@ OversetFE::solve() if (m_stage == 3) { // Compute diagnostics, e.g., residuals diag_computed = m_diag.compute( *d, m_u, m_un, m_bnorm, - m_symbcnodes, m_farfieldbcnodes ); + m_symbcnodes, m_farfieldbcnodes, + m_slipwallbcnodes ); // Increase number of iterations and physical time d->next(); // Advance physical time for local time stepping diff --git a/src/Inciter/OversetFE.hpp b/src/Inciter/OversetFE.hpp index dbb85d2177..6594d4acb3 100644 --- a/src/Inciter/OversetFE.hpp +++ b/src/Inciter/OversetFE.hpp @@ -222,7 +222,9 @@ class OversetFE : public CBase_OversetFE { p | m_bnormc; p | m_symbcnodes; p | m_farfieldbcnodes; + p | m_slipwallbcnodes; p | m_symbctri; + p | m_slipwallbctri; p | m_timedepbcnodes; p | m_timedepbcFn; p | m_stage; @@ -242,6 +244,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 @@ -335,8 +340,12 @@ class OversetFE : public CBase_OversetFE { std::unordered_set< std::size_t > m_symbcnodes; //! Unique set of nodes at which farfield BCs are set std::unordered_set< std::size_t > m_farfieldbcnodes; + //! Unique set of nodes at which slip wall BCs are set + std::unordered_set< std::size_t > m_slipwallbcnodes; //! Vector with 1 at symmetry BC boundary triangles std::vector< int > m_symbctri; + //! Vector with 1 at slip wall BC boundary triangles + std::vector< int > m_slipwallbctri; //! \brief Unique set of nodes at which time dependent BCs are set // for each time dependent BC std::vector< std::unordered_set< std::size_t > > m_timedepbcnodes; @@ -381,6 +390,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 +468,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:: diff --git a/src/PDE/CGPDE.hpp b/src/PDE/CGPDE.hpp index afbb6f6544..24c18fdfcc 100644 --- a/src/PDE/CGPDE.hpp +++ b/src/PDE/CGPDE.hpp @@ -162,6 +162,7 @@ class CGPDE { const std::pair< std::vector< std::size_t >, std::vector< std::size_t > >& esup, const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const std::vector< real >& vol, const std::vector< std::size_t >& edgenode, const std::vector< std::size_t >& edgeid, @@ -173,18 +174,18 @@ class CGPDE { real V, tk::Fields& R ) const { self->rhs( t, coord, inpoel, triinpoel, gid, bid, lid, dfn, psup, - esup, symbctri, vol, edgenode, edgeid, + esup, symbctri, slipwallbctri, vol, edgenode, edgeid, boxnodes, G, U, W, tp, V, R ); } //! Public interface to compute boundary surface integrals of pressure void bndPressureInt( const std::array< std::vector< real >, 3 >& coord, const std::vector< std::size_t >& triinpoel, - const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const tk::Fields& U, const std::array< tk::real, 3 >& CM, std::vector< real >& F ) const - { self->bndPressureInt( coord, triinpoel, symbctri, U, CM, F ); } + { self->bndPressureInt( coord, triinpoel, slipwallbctri, U, CM, F ); } //! Public interface for computing the minimum time step size real dt( const std::array< std::vector< real >, 3 >& coord, @@ -235,6 +236,17 @@ class CGPDE { const std::unordered_set< std::size_t >& nodes ) const { self->farfieldbc( U, coord, bnorm, nodes ); } + //! Public interface to set slip wall boundary conditions at nodes + void + slipwallbc( 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->slipwallbc( U, W, coord, bnorm, nodes ); } + //! Public interface to applying time dependent boundary conditions at nodes void timedepbc( tk::real t, @@ -347,6 +359,7 @@ class CGPDE { const std::pair< std::vector< std::size_t >, std::vector< std::size_t > >&, const std::vector< int >&, + const std::vector< int >&, const std::vector< real >&, const std::vector< std::size_t >&, const std::vector< std::size_t >&, @@ -397,6 +410,14 @@ class CGPDE { std::unordered_map< std::size_t, std::array< real, 4 > > >&, const std::unordered_set< std::size_t >& ) const = 0; + virtual void slipwallbc( + 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 > > >&, + const std::unordered_set< std::size_t >& ) const = 0; virtual void timedepbc( tk::real, tk::Fields&, @@ -474,6 +495,7 @@ class CGPDE { const std::pair< std::vector< std::size_t >, std::vector< std::size_t > >& esup, const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const std::vector< real >& vol, const std::vector< std::size_t >& edgenode, const std::vector< std::size_t >& edgeid, @@ -485,16 +507,16 @@ class CGPDE { real V, tk::Fields& R ) const override { data.rhs( t, coord, inpoel, triinpoel, gid, bid, lid, dfn, psup, - esup, symbctri, vol, edgenode, + esup, symbctri, slipwallbctri, vol, edgenode, edgeid, boxnodes, G, U, W, tp, V, R ); } void bndPressureInt( const std::array< std::vector< real >, 3 >& coord, const std::vector< std::size_t >& triinpoel, - const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const tk::Fields& U, const std::array< tk::real, 3 >& CM, std::vector< real >& F ) const override - { data.bndPressureInt( coord, triinpoel, symbctri, U, CM, F ); } + { data.bndPressureInt( coord, triinpoel, slipwallbctri, U, CM, F ); } real dt( const std::array< std::vector< real >, 3 >& coord, const std::vector< std::size_t >& inpoel, tk::real t, @@ -534,6 +556,15 @@ class CGPDE { std::array< real, 4 > > >& bnorm, const std::unordered_set< std::size_t >& nodes ) const override { data.farfieldbc( U, coord, bnorm, nodes ); } + void slipwallbc( + 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.slipwallbc( U, W, coord, bnorm, nodes ); } void timedepbc( tk::real t, diff --git a/src/PDE/CompFlow/CGCompFlow.hpp b/src/PDE/CompFlow/CGCompFlow.hpp index ae944f19d6..2fac287b56 100644 --- a/src/PDE/CompFlow/CGCompFlow.hpp +++ b/src/PDE/CompFlow/CGCompFlow.hpp @@ -402,6 +402,7 @@ class CompFlow { //! \param[in] psup Points surrounding points //! \param[in] esup Elements surrounding points //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles + //! \param[in] slipwallbctri Vector with 1 at slip BC boundary triangles //! \param[in] vol Nodal volumes //! \param[in] edgenode Local node IDs of edges //! \param[in] edgeid Edge ids in the order of access @@ -425,6 +426,7 @@ class CompFlow { const std::pair< std::vector< std::size_t >, std::vector< std::size_t > >& esup, const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const std::vector< real >& vol, const std::vector< std::size_t >& edgenode, const std::vector< std::size_t >& edgeid, @@ -455,7 +457,7 @@ class CompFlow { domainint( coord, gid, edgenode, edgeid, psup, dfn, U, W, Grad, R ); // compute boundary integrals - bndint( coord, triinpoel, symbctri, U, W, R ); + bndint( coord, triinpoel, symbctri, slipwallbctri, U, W, R ); // compute external (energy) sources const auto& icbox = g_inputdeck.get< tag::ic, tag::box >(); @@ -483,14 +485,14 @@ class CompFlow { //! Compute boundary pressure integrals (force) for rigid body motion //! \param[in] coord Mesh node coordinates //! \param[in] triinpoel Boundary triangle face connecitivity with local ids - //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles + //! \param[in] slipwallbctri Vector with 1 at symmetry BC boundary triangles //! \param[in] U Solution vector at recent time step //! \param[in] CM Center of mass //! \param[in,out] F Force vector (appended with torque vector) computed void bndPressureInt( const std::array< std::vector< real >, 3 >& coord, const std::vector< std::size_t >& triinpoel, - const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const tk::Fields& U, const std::array< tk::real, 3 >& CM, std::vector< real >& F ) const @@ -503,7 +505,7 @@ class CompFlow { // boundary integrals: compute surface integral of pressure (=force) for (std::size_t e=0; e, 3 >&, + 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 + { + // collect sidesets across all meshes + std::vector< std::size_t > swbc; + for (const auto& ibc : g_inputdeck.get< tag::bc >()) { + swbc.insert(swbc.end(), ibc.get< tag::slipwall >().begin(), + ibc.get< tag::slipwall >().end()); + } + + if (swbc.size() > 0) { // use slip bcs for this system + for (auto p : nodes) { // for all slipbc nodes + // for all user-def slipbc sets + for (std::size_t s=0; s(swbc[s])); + if (j != end(bnorm)) { + auto i = j->second.find(p); // find normal for node + if (i != end(j->second)) { + std::array< real, 3 > + n{ i->second[0], i->second[1], i->second[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 ); + // slip wall bc: 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]; + } + } + } + } + } + } + //! Apply user defined time dependent BCs //! \param[in] t Physical time //! \param[in,out] U Solution vector at recent time step @@ -1266,12 +1313,14 @@ class CompFlow { //! \param[in] coord Mesh node coordinates //! \param[in] triinpoel Boundary triangle face connecitivity with local ids //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles + //! \param[in] slipwallbctri Vector with 1 at slip wall BC boundary triangles //! \param[in] U Solution vector at recent time step //! \param[in] W Mesh velocity //! \param[in,out] R Right-hand side vector computed void bndint( const std::array< std::vector< real >, 3 >& coord, const std::vector< std::size_t >& triinpoel, const std::vector< int >& symbctri, + const std::vector< int >& slipwallbctri, const tk::Fields& U, const tk::Fields& W, tk::Fields& R ) const @@ -1325,30 +1374,34 @@ class CompFlow { real f[m_ncomp][3]; real p, vn; int sym = symbctri[e]; + int slip = slipwallbctri[e]; 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)); + vn = (sym || slip) ? 0.0 : (nx*(ruA/rA-w1A) + ny*(rvA/rA-w2A) + nz*(rwA/rA-w3A)); f[0][0] = rA*vn; f[1][0] = ruA*vn + p*nx; f[2][0] = rvA*vn + p*ny; f[3][0] = rwA*vn + p*nz; - f[4][0] = reA*vn + p*(sym ? 0.0 : (nx*ruA + ny*rvA + nz*rwA)/rA); + f[4][0] = reA*vn + p*((sym || slip) + ? 0.0 : (nx*ruA + ny*rvA + nz*rwA)/rA); p = m_mat_blk[0].compute< EOS::pressure >( rB, ruB/rB, rvB/rB, rwB/rB, reB ); - vn = sym ? 0.0 : (nx*(ruB/rB-w1B) + ny*(rvB/rB-w2B) + nz*(rwB/rB-w3B)); + vn = (sym || slip) ? 0.0 : (nx*(ruB/rB-w1B) + ny*(rvB/rB-w2B) + nz*(rwB/rB-w3B)); f[0][1] = rB*vn; f[1][1] = ruB*vn + p*nx; f[2][1] = rvB*vn + p*ny; f[3][1] = rwB*vn + p*nz; - f[4][1] = reB*vn + p*(sym ? 0.0 : (nx*ruB + ny*rvB + nz*rwB)/rB); + f[4][1] = reB*vn + p*((sym || slip) + ? 0.0 : (nx*ruB + ny*rvB + nz*rwB)/rB); p = m_mat_blk[0].compute< EOS::pressure >( rC, ruC/rC, rvC/rC, rwC/rC, reC ); - vn = sym ? 0.0 : (nx*(ruC/rC-w1C) + ny*(rvC/rC-w2C) + nz*(rwC/rC-w3C)); + vn = (sym || slip) ? 0.0 : (nx*(ruC/rC-w1C) + ny*(rvC/rC-w2C) + nz*(rwC/rC-w3C)); f[0][2] = rC*vn; f[1][2] = ruC*vn + p*nx; f[2][2] = rvC*vn + p*ny; f[3][2] = rwC*vn + p*nz; - f[4][2] = reC*vn + p*(sym ? 0.0 : (nx*ruC + ny*rvC + nz*rwC)/rC); + f[4][2] = reC*vn + p*((sym || slip) + ? 0.0 : (nx*ruC + ny*rvC + nz*rwC)/rC); // compute face area auto A6 = tk::area( x[N[0]], x[N[1]], x[N[2]], y[N[0]], y[N[1]], y[N[2]], diff --git a/src/PDE/CompFlow/DGCompFlow.hpp b/src/PDE/CompFlow/DGCompFlow.hpp index 0c46cfe76d..388a143504 100644 --- a/src/PDE/CompFlow/DGCompFlow.hpp +++ b/src/PDE/CompFlow/DGCompFlow.hpp @@ -78,13 +78,15 @@ class CompFlow { , invalidBC // Outlet BC not implemented , farfield , extrapolate - , invalidBC }, // No slip wall BC not implemented + , invalidBC // No slip wall BC not implemented + , symmetry }, // Slip equivalent to symmetry without mesh motion // BC Gradient functions { noOpGrad , noOpGrad , noOpGrad , noOpGrad , noOpGrad + , noOpGrad , noOpGrad } ) ); diff --git a/src/PDE/ConfigureCompFlow.cpp b/src/PDE/ConfigureCompFlow.cpp index 4992cb3c63..3804a681f7 100644 --- a/src/PDE/ConfigureCompFlow.cpp +++ b/src/PDE/ConfigureCompFlow.cpp @@ -142,6 +142,10 @@ infoCompFlow( std::map< ctr::PDEType, tk::ncomp_t >& cnt ) if (!sym.empty()) nfo.emplace_back( "Symmetry BC sideset(s)", parameters( sym ) ); + const auto& slip = ib.get< tag::slipwall >(); + if (!slip.empty()) + nfo.emplace_back( "Slip wall BC sideset(s)", parameters( slip ) ); + const auto& dir = ib.get< tag::dirichlet >(); if (!dir.empty()) nfo.emplace_back( "Dirichlet BC sideset(s)", parameters( dir ) ); diff --git a/src/PDE/MultiMat/DGMultiMat.hpp b/src/PDE/MultiMat/DGMultiMat.hpp index 908a713f04..a976bb3826 100644 --- a/src/PDE/MultiMat/DGMultiMat.hpp +++ b/src/PDE/MultiMat/DGMultiMat.hpp @@ -82,14 +82,16 @@ class MultiMat { , invalidBC // Outlet BC not implemented , farfield , extrapolate - , noslipwall }, + , noslipwall + , symmetry }, // Slip equivalent to symmetry without mesh motion // BC Gradient functions { noOpGrad , symmetryGrad , noOpGrad , noOpGrad , noOpGrad - , noOpGrad } + , noOpGrad + , symmetryGrad } ) ); // Inlet BC has a different structure than above BCs, so it must be diff --git a/src/PDE/MultiMat/FVMultiMat.hpp b/src/PDE/MultiMat/FVMultiMat.hpp index 9f1b112503..d04ee3338a 100644 --- a/src/PDE/MultiMat/FVMultiMat.hpp +++ b/src/PDE/MultiMat/FVMultiMat.hpp @@ -78,14 +78,16 @@ class MultiMat { , invalidBC // Outlet BC not implemented , farfield , extrapolate - , noslipwall }, + , noslipwall + , symmetry }, // Slip equivalent to symmetry without mesh motion // BC Gradient functions { noOpGrad , symmetryGrad , noOpGrad , noOpGrad , noOpGrad - , noOpGrad } + , noOpGrad + , symmetryGrad } ) ); // Inlet BC has a different structure than above BCs, so it must be diff --git a/src/PDE/MultiSpecies/DGMultiSpecies.hpp b/src/PDE/MultiSpecies/DGMultiSpecies.hpp index 0cc95fd9f1..4db379d307 100644 --- a/src/PDE/MultiSpecies/DGMultiSpecies.hpp +++ b/src/PDE/MultiSpecies/DGMultiSpecies.hpp @@ -81,14 +81,16 @@ class MultiSpecies { , invalidBC // Outlet BC not implemented , farfield , extrapolate - , noslipwall }, + , noslipwall + , symmetry }, // Slip equivalent to symmetry without mesh motion // BC Gradient functions { noOpGrad , symmetryGrad , noOpGrad , noOpGrad , noOpGrad - , noOpGrad } + , noOpGrad + , symmetryGrad } ) ); // EoS initialization diff --git a/src/PDE/Transport/CGTransport.hpp b/src/PDE/Transport/CGTransport.hpp index dd6bb33c62..62d5b7f932 100644 --- a/src/PDE/Transport/CGTransport.hpp +++ b/src/PDE/Transport/CGTransport.hpp @@ -227,6 +227,7 @@ class Transport { const std::pair< std::vector< std::size_t >, std::vector< std::size_t > >& esup, const std::vector< int >& symbctri, + const std::vector< int >&, const std::vector< real >& vol, const std::vector< std::size_t >&, const std::vector< std::size_t >& edgeid, @@ -413,6 +414,17 @@ class Transport { std::array< real, 4 > > >&, const std::unordered_set< std::size_t >& ) const {} + //! Set slip wall boundary conditions at nodes + void + slipwallbc( + tk::Fields&, + const tk::Fields&, + const std::array< std::vector< real >, 3 >&, + const std::unordered_map< int, + std::unordered_map< std::size_t, + std::array< real, 4 > > >&, + const std::unordered_set< std::size_t >& ) const {} + //! Apply user defined time dependent BCs (no-op for transport) void timedepbc( tk::real, diff --git a/src/PDE/Transport/DGTransport.hpp b/src/PDE/Transport/DGTransport.hpp index 4b4e0f0dfa..6c29c9142b 100644 --- a/src/PDE/Transport/DGTransport.hpp +++ b/src/PDE/Transport/DGTransport.hpp @@ -74,6 +74,7 @@ class Transport { , outlet , invalidBC // Characteristic BC not implemented , extrapolate + , invalidBC // Slip wall BC not implemented , invalidBC }, // No slip wall BC not implemented // BC Gradient functions { noOpGrad @@ -81,6 +82,7 @@ class Transport { , noOpGrad , noOpGrad , noOpGrad + , noOpGrad , noOpGrad } ) );