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, e.g., DiagCG.
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 : 65623 : match( [[maybe_unused]] tk::ctr::ncomp_t ncomp,
30 : : tk::real t,
31 : : tk::real dt,
32 : : const std::vector< tk::real >& tp,
33 : : const std::vector< tk::real >& dtp,
34 : : const tk::UnsMesh::Coords& coord,
35 : : const std::unordered_map< std::size_t, std::size_t >& lid,
36 : : const std::map< int, std::vector< std::size_t > >& bnode,
37 : : bool increment )
38 : : // *****************************************************************************
39 : : // Match user-specified boundary conditions at nodes for side sets
40 : : //! \param[in] ncomp Number of scalar components in PDE system
41 : : //! \param[in] t Physical time at which to query boundary conditions
42 : : //! \param[in] dt Time step size (for querying BC increments in time)
43 : : //! \param[in] tp Physical time for each mesh node
44 : : //! \param[in] dtp Time step size for each mesh node
45 : : //! \param[in] coord Mesh node coordinates
46 : : //! \param[in] lid Local node IDs associated to local node IDs
47 : : //! \param[in] bnode Map storing global mesh node IDs mapped to side set ids
48 : : //! \param[in] increment If true, evaluate the solution increment between
49 : : //! t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
50 : : //! \return Vector of pairs of bool and boundary condition value associated to
51 : : //! local mesh node IDs at which the user has set Dirichlet boundary
52 : : //! conditions for all systems of PDEs integrated. The bool indicates whether
53 : : //! the BC is set at the node for that component: if true, the real value is
54 : : //! the increment (from t to dt) in (or the value of) the BC specified for a
55 : : //! component.
56 : : //! \details Boundary conditions (BC), mathematically speaking, are applied on
57 : : //! finite surfaces. These finite surfaces are given by element sets (i.e., a
58 : : //! list of elements). This function queries Dirichlet boundary condition
59 : : //! values from all PDEs in the multiple systems of PDEs integrated at the
60 : : //! node lists associated to side set IDs, given by bnode. Each
61 : : //! PDE system returns a BC data structure. Note that the BC mesh nodes that
62 : : //! this function results in (stored in dirbc) only contains those nodes that
63 : : //! are supplied via bnode. i.e., in parallel only a part of the mesh is
64 : : //! worked on.
65 : : // *****************************************************************************
66 : : {
67 : : using inciter::g_cgpde;
68 : :
69 : : // Vector of pairs of bool and boundary condition value associated to mesh
70 : : // node IDs at which the user has set Dirichlet BCs for all PDEs integrated.
71 : : std::unordered_map< std::size_t,
72 : 65623 : std::vector< std::pair< bool, tk::real > > > dirbc;
73 : :
74 : : // Details for the algorithm below: PDE::dirbc() returns a new map that
75 : : // associates a vector of pairs associated to local node IDs. (The pair is a
76 : : // pair of bool and real value, the former is the fact that the BC is to be
77 : : // set while the latter is the value if it is to be set). The length of this
78 : : // NodeBC vector, returning from each system of PDEs equals to the number of
79 : : // scalar components the given PDE integrates. Here we contatenate this map
80 : : // for all PDEs being integrated. If there are multiple BCs set at a mesh node
81 : : // (dirbc::key), either because (1) in the same PDE system the user prescribed
82 : : // BCs on side sets that share nodes or (2) because more than a single PDE
83 : : // system assigns BCs to a given node (on different variables), the NodeBC
84 : : // vector must be correctly stored. "Correctly" here means that the size of
85 : : // the NodeBC vectors must all be the same and equal to the sum of all scalar
86 : : // components integrated by all PDE systems integrated. Example: single-phase
87 : : // compressible flow (density, momentum, energy = 5) + transported scalars of
88 : : // 10 variables -> NodeBC vector length = 15. Note that in case (1) above a
89 : : // new node encountered must "overwrite" the already existing space for the
90 : : // NodeBC vector. "Overwrite" here means that it should keep the existing BCs
91 : : // and add the new ones yielding the union the two prescription for BCs but in
92 : : // the same space that already exist in the NodeBC vector. In case (2),
93 : : // however, the NodeBC pairs must go to the location in the vector assigned to
94 : : // the given PDE system, i.e., using the above example BCs for the 10 (or
95 : : // less) scalars should go in the positions starting at 5, leaving the first 5
96 : : // false, indicating no BCs for the flow variables.
97 : : //
98 : : // Note that the logic described above is only partially implemented at this
99 : : // point. What works is the correct insertion of multiple BCs for nodes shared
100 : : // among multiple side sets, e.g., corners, originating from the same PDE
101 : : // system. What is not yet implemented is the case when there are no BCs set
102 : : // for flow variables but there are BCs for transport, the else branch below
103 : : // will incorrectly NOT skip the space for the flow variables. In other words,
104 : : // this only works for a single PDE system and a sytem of systems. This
105 : : // machinery is only tested with a single system of PDEs at this point.
106 : : //
107 : : // When a particular node belongs to two or more side sets with different BCs,
108 : : // there is an ambiguity as to which of the multiple BCs should be applied to
109 : : // the node. This issue is described in case (1) above. In the current
110 : : // implementation, every side set applies the BC to the common node in
111 : : // question, successively overwriting the BC applied by the previous side set.
112 : : // Effectively, the BC corresponding to the last side set ID is applied to the
113 : : // common node. Since bnode is an ordered map, the side set with a larger
114 : : // id wins if a node belongs to multiple side sets.
115 : :
116 : : // Lambda to convert global to local node ids of a list of nodes
117 : 10989591 : auto local = [ &lid ]( const std::vector< std::size_t >& gnodes ){
118 [ + - ]: 111917 : std::vector< std::size_t > lnodes( gnodes.size() );
119 [ + + ]: 10989591 : for (std::size_t i=0; i<gnodes.size(); ++i)
120 [ + - ]: 10877674 : lnodes[i] = tk::cref_find( lid, gnodes[i] );
121 : 111917 : return lnodes;
122 : 65623 : };
123 : :
124 : : // Query Dirichlet BCs for all PDEs integrated and assign to nodes
125 [ + + ]: 177540 : for (const auto& s : bnode) { // for all side sets passed in
126 : 111917 : std::size_t c = 0;
127 [ + - ]: 223834 : auto l = local(s.second); // generate local node ids on side set
128 [ + + ]: 223834 : for (std::size_t eq=0; eq<g_cgpde.size(); ++eq) {
129 : : // query Dirichlet BCs at nodes of this side set
130 : : auto eqbc =
131 [ + - ][ + - ]: 223834 : g_cgpde[eq].dirbc( t, dt, tp, dtp, {s.first,l}, coord, increment );
132 [ + + ]: 1473191 : for (const auto& n : eqbc) {
133 : 1361274 : auto id = n.first; // BC node ID
134 : 1361274 : const auto& bcs = n.second; // BCs
135 [ + - ]: 1361274 : auto& nodebc = dirbc[ id ]; // BCs to be set for node
136 [ + + ]: 1361274 : if (nodebc.size() > c) { // node already has BCs from this PDE
137 [ - + ][ - - ]: 267006 : Assert( nodebc.size() == c+bcs.size(), "Size mismatch" );
[ - - ][ - - ]
138 [ + + ]: 1545836 : for (std::size_t i=0; i<bcs.size(); i++) {
139 [ + - ]: 1278830 : 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 [ + - ]: 1094268 : nodebc.insert( end(nodebc), begin(bcs), end(bcs) );
145 : : }
146 : : }
147 [ + + ]: 111917 : if (!eqbc.empty()) c += eqbc.cbegin()->second.size();
148 : : }
149 : : }
150 : :
151 : : // Verify the size of each NodeBC vectors. They must have the same lengths and
152 : : // equal to the total number of scalar components for all systems of PDEs
153 : : // integrated.
154 [ + - ][ - + ]: 1159891 : Assert( std::all_of( begin(dirbc), end(dirbc),
[ - - ][ - - ]
[ - - ]
155 : : [ ncomp ]( const auto& n ){ return n.second.size() == ncomp; } ),
156 : : "Size of NodeBC vector incorrect" );
157 : :
158 : 131246 : return dirbc;
159 : : }
160 : :
161 : : bool
162 : 18283 : correctBC( const tk::Fields& a,
163 : : const tk::Fields& dul,
164 : : const std::unordered_map< std::size_t,
165 : : std::vector< std::pair< bool, tk::real > > >& dirbc )
166 : : // *****************************************************************************
167 : : // Verify that the change in the solution at those nodes where Dirichlet
168 : : // boundary conditions are set is exactly the amount the BCs prescribe
169 : : //! \param[in] a Limited antidiffusive element contributions (from FCT)
170 : : //! \param[in] dul Low order solution increment
171 : : //! \param[in] dirbc Vector of boundary conditions (true/false + BC value) for
172 : : //! all scalar components integrated associated of all systems to local node
173 : : //! ID
174 : : //! \return True if solution is correct at Dirichlet boundary condition nodes
175 : : //! \details We loop through the map that associates a vector of of boundary
176 : : //! conditions (true/false, indicating whether the BC is set + BC value if
177 : : //! true) for all scalar components integrated associated of all systems to
178 : : //! global node IDs. Then for all scalar components of all systems of systems
179 : : //! of PDEs integrated if a BC is to be set for a given component, we compute
180 : : //! the low order solution increment + the anti-diffusive element
181 : : //! contributions (in FCT), which is the current solution increment (to be
182 : : //! used to update the solution at time n in FCT) at that node. This solution
183 : : //! increment must equal the BC prescribed at the given node as we solve for
184 : : //! solution increments. If not, the BCs are not set correctly, which is an
185 : : //! error.
186 : : // *****************************************************************************
187 : : {
188 [ + + ]: 528007 : for (const auto& [i,bc] : dirbc) {
189 [ - + ][ - - ]: 509724 : Assert( bc.size() == dul.nprop(), "Size mismatch" );
[ - - ][ - - ]
190 [ + + ]: 2558044 : for (std::size_t c=0; c<bc.size(); ++c) {
191 [ + - ][ - + ]: 4096640 : if ( bc[c].first &&
[ - + ]
192 [ + - ][ + - ]: 2048320 : std::abs( dul(i,c,0) + a(i,c,0) - bc[c].second ) >
193 : 2048320 : std::numeric_limits< tk::real >::epsilon() )
194 : : {
195 : 0 : return false;
196 : : }
197 : : }
198 : : }
199 : :
200 : 18283 : return true;
201 : : }
202 : :
203 : : } // inciter::
|
