Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Inciter/NodeBC.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Boundary conditions for nodal discretizations
       9                 :            :   \details   Boundary conditions for nodal discretizations, such as continuous
      10                 :            :     Galerkin finite elements.
      11                 :            : */
      12                 :            : // *****************************************************************************
      13                 :            : 
      14                 :            : #include <vector>
      15                 :            : #include <map>
      16                 :            : #include <unordered_map>
      17                 :            : #include <algorithm>
      18                 :            : 
      19                 :            : #include "NodeBC.hpp"
      20                 :            : #include "CGPDE.hpp"
      21                 :            : #include "Fields.hpp"
      22                 :            : #include "Vector.hpp"
      23                 :            : 
      24                 :            : namespace inciter {
      25                 :            : 
      26                 :            : extern std::vector< CGPDE > g_cgpde;
      27                 :            : 
      28                 :            : std::unordered_map< std::size_t, std::vector< std::pair< bool, tk::real > > >
      29                 :      45114 : match( std::size_t meshid,
      30                 :            :        [[maybe_unused]] tk::ncomp_t ncomp,
      31                 :            :        tk::real t,
      32                 :            :        tk::real dt,
      33                 :            :        const std::vector< tk::real >& tp,
      34                 :            :        const std::vector< tk::real >& dtp,
      35                 :            :        const tk::UnsMesh::Coords& coord,
      36                 :            :        const std::unordered_map< std::size_t, std::size_t >& lid,
      37                 :            :        const std::map< int, std::vector< std::size_t > >& bnode,
      38                 :            :        bool increment )
      39                 :            : // *****************************************************************************
      40                 :            : //  Match user-specified boundary conditions at nodes for side sets
      41                 :            : //! \param[in] ncomp Number of scalar components in PDE system
      42                 :            : //! \param[in] t Physical time at which to query boundary conditions
      43                 :            : //! \param[in] dt Time step size (for querying BC increments in time)
      44                 :            : //! \param[in] tp Physical time for each mesh node
      45                 :            : //! \param[in] dtp Time step size for each mesh node
      46                 :            : //! \param[in] coord Mesh node coordinates
      47                 :            : //! \param[in] lid Local node IDs associated to local node IDs
      48                 :            : //! \param[in] bnode Map storing global mesh node IDs mapped to side set ids
      49                 :            : //! \param[in] increment If true, evaluate the solution increment between
      50                 :            : //!   t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
      51                 :            : //! \return Vector of pairs of bool and boundary condition value associated to
      52                 :            : //!   local mesh node IDs at which the user has set Dirichlet boundary
      53                 :            : //!   conditions for all systems of PDEs integrated. The bool indicates whether
      54                 :            : //!   the BC is set at the node for that component: if true, the real value is
      55                 :            : //!   the increment (from t to dt) in (or the value of) the BC specified for a
      56                 :            : //!   component.
      57                 :            : //! \details Boundary conditions (BC), mathematically speaking, are applied on
      58                 :            : //!   finite surfaces. These finite surfaces are given by element sets (i.e., a
      59                 :            : //!   list of elements). This function queries Dirichlet boundary condition
      60                 :            : //!   values from all PDEs in the multiple systems of PDEs integrated at the
      61                 :            : //!   node lists associated to side set IDs, given by bnode. Each
      62                 :            : //!   PDE system returns a BC data structure. Note that the BC mesh nodes that
      63                 :            : //!   this function results in (stored in dirbc) only contains those nodes that
      64                 :            : //!   are supplied via bnode. i.e., in parallel only a part of the mesh is
      65                 :            : //!   worked on.
      66                 :            : // *****************************************************************************
      67                 :            : {
      68                 :            :   using inciter::g_cgpde;
      69                 :            : 
      70                 :            :   // Vector of pairs of bool and boundary condition value associated to mesh
      71                 :            :   // node IDs at which the user has set Dirichlet BCs for all PDEs integrated.
      72                 :            :   std::unordered_map< std::size_t,
      73                 :            :     std::vector< std::pair< bool, tk::real > > > dirbc;
      74                 :            : 
      75                 :            :   // Details for the algorithm below: PDE::dirbc() returns a new map that
      76                 :            :   // associates a vector of pairs associated to local node IDs. (The pair is a
      77                 :            :   // pair of bool and real value, the former is the fact that the BC is to be
      78                 :            :   // set while the latter is the value if it is to be set). The length of this
      79                 :            :   // NodeBC vector, returning from each system of PDEs equals to the number of
      80                 :            :   // scalar components the given PDE integrates. Here we contatenate this map
      81                 :            :   // for all PDEs being integrated. If there are multiple BCs set at a mesh node
      82                 :            :   // (dirbc::key), either because (1) in the same PDE system the user prescribed
      83                 :            :   // BCs on side sets that share nodes or (2) because more than a single PDE
      84                 :            :   // system assigns BCs to a given node (on different variables), the NodeBC
      85                 :            :   // vector must be correctly stored. "Correctly" here means that the size of
      86                 :            :   // the NodeBC vectors must all be the same and equal to the sum of all scalar
      87                 :            :   // components integrated by all PDE systems integrated. Example: single-phase
      88                 :            :   // compressible flow (density, momentum, energy = 5) + transported scalars of
      89                 :            :   // 10 variables -> NodeBC vector length = 15. Note that in case (1) above a
      90                 :            :   // new node encountered must "overwrite" the already existing space for the
      91                 :            :   // NodeBC vector. "Overwrite" here means that it should keep the existing BCs
      92                 :            :   // and add the new ones yielding the union the two prescription for BCs but in
      93                 :            :   // the same space that already exist in the NodeBC vector. In case (2),
      94                 :            :   // however, the NodeBC pairs must go to the location in the vector assigned to
      95                 :            :   // the given PDE system, i.e., using the above example BCs for the 10 (or
      96                 :            :   // less) scalars should go in the positions starting at 5, leaving the first 5
      97                 :            :   // false, indicating no BCs for the flow variables.
      98                 :            :   //
      99                 :            :   // Note that the logic described above is only partially implemented at this
     100                 :            :   // point. What works is the correct insertion of multiple BCs for nodes shared
     101                 :            :   // among multiple side sets, e.g., corners, originating from the same PDE
     102                 :            :   // system. What is not yet implemented is the case when there are no BCs set
     103                 :            :   // for flow variables but there are BCs for transport, the else branch below
     104                 :            :   // will incorrectly NOT skip the space for the flow variables. In other words,
     105                 :            :   // this only works for a single PDE system and a sytem of systems. This
     106                 :            :   // machinery is only tested with a single system of PDEs at this point.
     107                 :            :   //
     108                 :            :   // When a particular node belongs to two or more side sets with different BCs,
     109                 :            :   // there is an ambiguity as to which of the multiple BCs should be applied to
     110                 :            :   // the node. This issue is described in case (1) above. In the current
     111                 :            :   // implementation, every side set applies the BC to the common node in
     112                 :            :   // question, successively overwriting the BC applied by the previous side set.
     113                 :            :   // Effectively, the BC corresponding to the last side set ID is applied to the
     114                 :            :   // common node. Since bnode is an ordered map, the side set with a larger
     115                 :            :   // id wins if a node belongs to multiple side sets.
     116                 :            : 
     117                 :            :   // Lambda to convert global to local node ids of a list of nodes
     118                 :    8981959 :   auto local = [ &lid ]( const std::vector< std::size_t >& gnodes ){
     119                 :      73684 :     std::vector< std::size_t > lnodes( gnodes.size() );
     120         [ +  + ]:    9055643 :     for (std::size_t i=0; i<gnodes.size(); ++i)
     121                 :    8981959 :       lnodes[i] = tk::cref_find( lid, gnodes[i] );
     122                 :      73684 :     return lnodes;
     123                 :      45114 :   };
     124                 :            : 
     125                 :            :   // Query Dirichlet BCs for all PDEs integrated and assign to nodes
     126         [ +  + ]:     118798 :   for (const auto& s : bnode) {     // for all side sets passed in
     127                 :            :     std::size_t c = 0;
     128         [ +  - ]:      73684 :     auto l = local(s.second);   // generate local node ids on side set
     129                 :            :     // query Dirichlet BCs at nodes of this side set
     130                 :            :     auto eqbc =
     131 [ +  - ][ -  - ]:     147368 :       g_cgpde[meshid].dirbc( t, dt, tp, dtp, {s.first,l}, coord, increment );
     132         [ +  + ]:     856881 :     for (const auto& n : eqbc) {
     133         [ +  - ]:     783197 :       auto id = n.first;                      // BC node ID
     134                 :            :       const auto& bcs = n.second;             // BCs
     135                 :            :       auto& nodebc = dirbc[ id ];     // BCs to be set for node
     136         [ +  + ]:     783197 :       if (nodebc.size() > c) {        // node already has BCs from this PDE
     137                 :            :         Assert( nodebc.size() == c+bcs.size(), "Size mismatch" );
     138         [ +  + ]:     923016 :         for (std::size_t i=0; i<bcs.size(); i++) {
     139         [ +  - ]:     769180 :           if (bcs[i].first) nodebc[c+i] = bcs[i];
     140                 :            :         }
     141                 :            :       } else {        // node does not yet have BCs from this PDE
     142                 :            :         // This branch needs to be completed for system of systems of PDEs.
     143                 :            :         // See note above.
     144         [ +  - ]:     629361 :         nodebc.insert( end(nodebc), begin(bcs), end(bcs) );
     145                 :            :       }
     146                 :            :     }
     147                 :            :     if (!eqbc.empty()) c += eqbc.cbegin()->second.size();
     148                 :            :   }
     149                 :            : 
     150                 :            :   // Verify the size of each NodeBC vectors. They must have the same lengths and
     151                 :            :   // equal to the total number of scalar components for all systems of PDEs
     152                 :            :   // integrated.
     153                 :            :   Assert( std::all_of( begin(dirbc), end(dirbc),
     154                 :            :             [ ncomp ]( const auto& n ){ return n.second.size() == ncomp; } ),
     155                 :            :           "Size of NodeBC vector incorrect" );
     156                 :            :  
     157                 :      45114 :   return dirbc;
     158                 :            : }
     159                 :            : 
     160                 :            : } // inciter::