Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Control/SystemComponents.hpp
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 Operations on numbers of scalar components of systems of equations
9 : : \details Operations on numbers of scalar components of systems of equations,
10 : : e.g. multiple equation sets of a physical model or a set of (stochastic
11 : : ordinary or partial differential) equations of different kinds.
12 : :
13 : : _Problem:_ We are given a type that is a tk::tuple::tagged_tuple that
14 : : contains an arbitrary number of std::vectors of integers. The number of
15 : : vectors are fixed at compile-time (accessed via tags) but their (differing)
16 : : length is only known at run-time (after parsing user input). What we need
17 : : are functions that operate on this data structure and return, e.g., the
18 : : total number of components in the whole system or the offset in the whole
19 : : data structure for a given tag. The functions should thus be able to operate
20 : : on a list of types, i.e., a double for-loop over all tags and associated
21 : : vectors - one at compile-time and the other one at run-time.
22 : :
23 : : _Example:_ An example is to define storage for systems of stochastic
24 : : differential equations as in walker::ctr::ncomps. In
25 : : [walker](walker.html) various types of stochastic differential
26 : : equations are implemented and numerically integrated in time in order to
27 : : compute and learn about their statistical behavior. Since different types of
28 : : systems can be parameterized differently, there can be several systems of
29 : : the same type (the number is user-defined so only known at run-time). For
30 : : example, it is possible to numerically integrate two systems of Dirichlet
31 : : equations, see DiffEq/Dirichlet.h, one with 4, the other one with 5 scalar
32 : : variables, parameterized differently, and estimate their arbitrary coupled
33 : : statistics and joint PDFs. This requires that each type of equation system
34 : : has a vector of integers storing the number of scalar variables.
35 : :
36 : : _Solution:_ Looping through a list of types is done using the [Brigand]
37 : : (https://github.com/edouarda/brigand)'s [MetaProgramming Library].
38 : : Such operations on types happen at compile-time, i.e., the code runs inside
39 : : the compiler and only its result gets compiled into code to be run at
40 : : run-time. Advantages are abstraction and generic code that is independent
41 : : of the size and order of the tags in the tuple (and the associated vectors).
42 : : Though it is not of primary concern for this data structure, this solution
43 : : also generates efficient code, since part of the algorithm (for computing
44 : : the offset, the total number of components, etc.) runs inside the compiler.
45 : : All of this is kept generic, so it does not need to be changed when a new
46 : : pair of tag and system of equations are added to the underlying tuple.
47 : : Thus, _the main advantage is that changing the order of the vectors in the
48 : : tuple or adding new ones at arbitrary locations will not require change to
49 : : client-code._
50 : :
51 : : _An alternative approach:_ Of course this could have been simply done with a
52 : : purely run-time vector of vector of integers. However, that would not have
53 : : allowed a tag-based access to the different systems of equations. (Sure, one
54 : : can define an enum for equation indexing, but that just adds code to what is
55 : : already need to be changed if a change happens to the underlying data
56 : : structure.) As a consequence, a component in a vector of vector would have
57 : : had to be accessed as something like, ncomps[2][3], which is more
58 : : error-prone to read than ncomps< tag::dirichlet >( 3 ). Furthermore, that
59 : : design would have resulted in double-loops at run-time, instead of letting
60 : : the compiler loop over the types at compile-time and do the run-time loops
61 : : only for what is only known at run-time. Additionally, and most importantly,
62 : : _any time a new system is introduced (or the order is changed), all
63 : : client-code would have to be changed._
64 : :
65 : : For example client-code, see walker::Dirichlet::Dirichlet() in
66 : : DiffEq/Dirichlet.h, or walker::Integrator::Integrator in Walker/Integrator.h.
67 : : */
68 : : // *****************************************************************************
69 : : #ifndef SystemComponents_h
70 : : #define SystemComponents_h
71 : :
72 : : #include <vector>
73 : : #include <algorithm>
74 : : #include <numeric>
75 : : #include <string>
76 : :
77 : : #include <brigand/sequences/list.hpp>
78 : : #include <brigand/algorithms/wrap.hpp>
79 : :
80 : : #include "NoWarning/flatten.hpp"
81 : : #include "NoWarning/transform.hpp"
82 : :
83 : : #include "TaggedTuple.hpp"
84 : : #include "StatCtr.hpp"
85 : : #include "Keywords.hpp"
86 : : #include "Tags.hpp"
87 : :
88 : : namespace tk {
89 : : //! Toolkit control, general purpose user input to internal data transfer
90 : : namespace ctr {
91 : :
92 : : //! Inherit type of number of components from keyword 'ncomp'
93 : : using ncomp_t = kw::ncomp::info::expect::type;
94 : :
95 : : //! \brief Map associating offsets to dependent variables for systems
96 : : //! \details This map associates offsets of systems of differential
97 : : //! equations in a larger data array storing dependent variables for all
98 : : //! scalar components of a system of systems. These offsets are where a
99 : : //! particular system starts and their field (or component) ids then can be
100 : : //! used to access an individual scalar component based on theses offsets.
101 : : //! \note We use a case-insensitive character comparison functor for the
102 : : //! offset map, since the keys (dependent variables) in the offset map are
103 : : //! only used to indicate the equation's dependent variable, however, queries
104 : : //! (using std::map's member function, find) can be fired up for both ordinary
105 : : //! and central moments of dependent variables (which are denoted by upper and
106 : : //! lower case, characters, respectively) for which the offset (for the same
107 : : //! dependent variable) should be the same.
108 : : using OffsetMap = std::map< char, ncomp_t, CaseInsensitiveCharLess >;
109 : :
110 : : //! \brief Map associating number of scalar components to dependent variables
111 : : //! for systems
112 : : //! \details This map associates the number of properties (scalar components)
113 : : //! of systems of differential equations for all scalar components of a
114 : : //! system of systems.
115 : : //! \note We use a case-insensitive character comparison functor to be
116 : : //! consistent with OffsetMap.
117 : : using NcompMap = std::map< char, ncomp_t, CaseInsensitiveCharLess >;
118 : :
119 : : //! Helper for converting a brigand::list to a tagged_tuple
120 : : template< typename... T >
121 : : using tagged_tuple_wrapper = typename tk::TaggedTuple< brigand::list<T...> >;
122 : :
123 : : //! Helper for converting a brigand::list to a tagged_tuple
124 : : template< typename L >
125 : : using as_tagged_tuple = brigand::wrap< L, tagged_tuple_wrapper >;
126 : :
127 : : //! Number of components storage as a vector for a system of equations
128 : : //! \details This is only a helper class, defining a type 'type' for
129 : : //! brigand::apply, so it can be used for defining a base for ncomponents
130 : : struct ComponentVector : public std::vector< ncomp_t > {
131 : : using type = std::vector< ncomp_t >;
132 : : };
133 : :
134 : : //! \brief Number of components storage
135 : : //! \details All this trickery with template meta-programming allows the code
136 : : //! below to be generic. As a result, adding a new component requires adding a
137 : : //! single line (a tag and its type) to the already existing list, see
138 : : //! typedefs 'ncomps'. The member functions, doing initialization, computing
139 : : //! the number of total components, the offset for a given tag, and computing
140 : : //! the offset map, need no change -- even if the order of the number of
141 : : //! components change.
142 : : template< typename... Tags >
143 : : class ncomponents : public
144 : : // tk::tuple::tagged_tuple< tag1, vec1, tag2, vec2, ... >
145 : : as_tagged_tuple< brigand::flatten< brigand::transform< brigand::list<Tags...>,
146 : : brigand::bind< brigand::list, brigand::_1, ComponentVector > > > > {
147 : :
148 : : private:
149 : : //! Access vector of number of components of an eq system as const-ref
150 : : template< typename Eq >
151 : 44091 : constexpr const auto& vec() const { return this->template get< Eq >(); }
152 : :
153 : : //! Access vector of number of components of an eq system as reference
154 : : template< typename Eq >
155 : 20988 : constexpr auto& vec() { return this->template get< Eq >(); }
156 : :
157 : : public:
158 : : //! Default constructor: set defaults to zero for all number of components
159 : 662 : ncomponents() {
160 [ + - ][ + - ]: 662 : ( ... , std::fill( begin(vec<Tags>()), end(vec<Tags>()), 0 ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
161 : 662 : }
162 : :
163 : : //! Compute total number of components in the systems of systems configured
164 : : //! \return Total number of components
165 : 957 : ncomp_t nprop() const noexcept {
166 : 957 : return (... + std::accumulate(begin(vec<Tags>()), end(vec<Tags>()), 0u));
167 : : }
168 : :
169 : : //! \brief Compute the offset for a given equation tag, i.e., the sum of the
170 : : //! number of components up to a given tag
171 : : //! \return offset for tag
172 : : //! \param[in] c Index for system given by template argument tag
173 : : //! \details This offset is used to index into the data array (the equation
174 : : //! systems operate on during the numerical solution) and get to the
175 : : //! beginning of data for a given differential equation system.
176 : : template< typename tag >
177 : 1638 : ncomp_t offset( ncomp_t c ) const noexcept {
178 : 1638 : ncomp_t offset = 0;
179 : 1638 : bool found = false;
180 : 27734 : ( ... , [&](){
181 : : if (std::is_same_v< tag, Tags >) {
182 : 1638 : const auto& v = vec<Tags>();
183 : : // Make sure we are not trying to index beyond the length for this eq
184 [ - + ][ - - ]: 1638 : Assert( v.size() >= c, "Indexing out of bounds" );
[ - - ][ - - ]
[ - + ][ - - ]
[ - - ][ - - ]
[ - + ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
185 : : // If we have found the tag we are looking for, we count up to the
186 : : // given system index and add those to the offset
187 [ + + ][ + + ]: 1642 : for (ncomp_t q=0; q<c; ++q) offset += v[q];
[ - + ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
188 : 1638 : found = true;
189 [ + + ][ + + ]: 24458 : } else if (!found) {
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ - + ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ - + ]
[ - + ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
190 : : // If we have not found the tag we are looking for, we add all the
191 : : // number of scalars for that tag to the offset
192 : 11857 : const auto& v = vec<Tags>();
193 : 11857 : offset = std::accumulate( begin(v), end(v), offset );
194 : 1638 : } }() );
195 : 1638 : return offset;
196 : : }
197 : :
198 : : //! Compute map of offsets associated to dependent variables
199 : : //! \param[in] d Input deck to operate on
200 : : //! \return Map of offsets associated to dependent variables
201 : : template< class InputDeck >
202 : 957 : OffsetMap offsetmap( const InputDeck& d ) const {
203 : 957 : OffsetMap map;
204 : 15298 : ( ... , [&](){
205 : 15298 : const auto& depvar = d.template get< tag::param, Tags, tag::depvar >();
206 : 15298 : ncomp_t c = 0;
207 : 15298 : const auto& ncomps = d.template get< tag::component >();
208 [ + + ][ + + ]: 16438 : for (auto v : depvar)
[ + + ][ + + ]
[ - + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ + + ][ + + ]
[ - + ][ + + ]
[ + + ][ + + ]
209 [ + - ][ + - ]: 2097 : map[ v ] = ncomps.template offset< Tags >( c++ ); }() );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ + - ]
[ + - ][ + - ]
210 : 957 : return map;
211 : : }
212 : :
213 : : //! \brief Compute map of number of properties (scalar components)
214 : : //! associated to dependent variables
215 : : //! \param[in] d Input deck to operate on
216 : : //! \return Map of number of properties associated to dependent variables
217 : : template< class InputDeck >
218 : 1 : NcompMap ncompmap( const InputDeck& d ) const {
219 : 1 : NcompMap map;
220 : 2 : ( ... , [&](){
221 : 2 : const auto& depvar = d.template get< tag::param, Tags, tag::depvar >();
222 : 2 : const auto& ncomps = d.template get< tag::component >();
223 : 2 : const auto& ncvec = ncomps.template get<Tags>();
224 [ - + ][ - - ]: 2 : Assert( ncvec.size() == depvar.size(), "ncompsize != depvarsize" );
[ - - ][ - - ]
[ - + ][ - - ]
[ - - ][ - - ]
225 : 2 : ncomp_t c = 0;
226 [ + - ][ + - ]: 7 : for (auto v : depvar) map[ v ] = ncvec[c++]; }() );
[ + + ][ + - ]
[ + + ][ + - ]
227 : 1 : return map;
228 : : }
229 : :
230 : : //! \brief Return vector of dependent variables + component id for all
231 : : //! equations configured
232 : : //! \param[in] deck Input deck to operate on
233 : : //! \return Vector of dependent variables + comopnent id for all equations
234 : : //! configured. The length of this vector equals the total number of
235 : : //! components configured, see nprop(), containing the depvar + the
236 : : //! component index relative to the given equation. E.g., c1, c2, u1, u2,
237 : : //! u3, u4, u5.
238 : : template< class InputDeck >
239 : : std::vector< std::string > depvar( const InputDeck& deck ) const {
240 : : std::vector< std::string > d;
241 : : ( ..., [&](){
242 : : const auto& dveq = deck.template get< tag::param, Tags, tag::depvar >();
243 : : const auto& nceq = deck.template get< tag::component, Tags >();
244 : : Assert( dveq.size() == nceq.size(), "Size mismatch" );
245 : : std::size_t e = 0;
246 : : for (auto v : dveq) {
247 : : for (std::size_t c=0; c<nceq[e]; ++c )
248 : : d.push_back( v + std::to_string(c+1) );
249 : : ++e;
250 : : } }() );
251 : : return d;
252 : : }
253 : : };
254 : :
255 : : } // ctr::
256 : : } // tk::
257 : :
258 : : #endif // SystemComponents_h
|
