Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Control/StatCtr.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Types and associated functions to deal with moments and PDFs
10 : : \details Types and associated functions to deal with statistical moments and
11 : : probability density functions.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef StatControl_h
15 : : #define StatControl_h
16 : :
17 : : #include "Types.hpp"
18 : : #include "Exception.hpp"
19 : : #include "PUPUtil.hpp"
20 : :
21 : : namespace tk {
22 : : namespace ctr {
23 : :
24 : : //! \brief Moment specifies which type of moment is computed for a quantity in
25 : : //! a Term
26 : : enum class Moment : uint8_t { ORDINARY=0, //!< Full variable
27 : : CENTRAL //!< Fluctuation
28 : : };
29 : :
30 : : //! \brief Term is a Moment of a quantity with a field ID to be ensemble
31 : : //! averaged
32 : : //! \details Internally the numbering of field IDs starts from 0, but presented
33 : : //! to the user, e.g., in screen-output, as starting from 1.
34 : : struct Term {
35 : : char var; //!< Variable name
36 : : uint64_t field; //!< Field ID
37 : : Moment moment; //!< Moment type: ordinary, central
38 : :
39 : : /** @name Pack/Unpack: Serialize Term object for Charm++ */
40 : : ///@{
41 : : //! Pack/Unpack serialize member function
42 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
43 : : void pup( PUP::er& p ) {
44 : : p | var;
45 : : p | field;
46 : : PUP::pup( p, moment );
47 : : }
48 : : //! \brief Pack/Unpack serialize operator|
49 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
50 : : //! \param[in,out] t Term object reference
51 : : friend void operator|( PUP::er& p, Term& t ) { t.pup(p); }
52 : : ///@}
53 : :
54 : : //! \brief Constructor: initialize all state data
55 : : //! \param[in] v Variable name
56 : : //! \param[in] f Field ID
57 : : //! \param[in] m Moment type enum: Moment::ORDINARY or Moment::CENTRAL
58 : 12 : explicit Term( char v = 0, uint64_t f = 0, Moment m = Moment::ORDINARY ) :
59 : 12 : var( v ), field( f ), moment( m ) {}
60 : :
61 : : //! \brief Equal operator for, e.g., finding unique elements, used by, e.g.,
62 : : //! std::unique().
63 : : //! \details Test on field, moment, and var
64 : : //! \param[in] term Term to compare
65 : : //! \return Boolean indicating if term equals 'this'
66 : 14 : bool operator== ( const Term& term ) const {
67 [ + - ][ + - ]: 14 : if (var == term.var && field == term.field && moment == term.moment)
[ + - ]
68 : 14 : return true;
69 : : else
70 : 0 : return false;
71 : : }
72 : :
73 : : //! \brief Less-than operator for ordering, used by, e.g., std::sort().
74 : : //! \details Test on var, field, and moment.
75 : : //! \param[in] term Term to compare
76 : : //! \return Boolean indicating if term is less than 'this'
77 : 32 : bool operator< ( const Term& term ) const {
78 [ + + ]: 32 : if (var < term.var)
79 : 8 : return true;
80 [ + + ][ - + ]: 24 : else if (var == term.var && field < term.field)
81 : 0 : return true;
82 [ + + ][ + - ]: 24 : else if (var == term.var && field == term.field && moment < term.moment)
[ - + ]
83 : 0 : return true;
84 : : else
85 : 24 : return false;
86 : : }
87 : : };
88 : :
89 : : //! \brief Pack/Unpack: Namespace-scope serialize Term object for Charm++
90 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
91 : : //! \param[in,out] t Term object reference
92 : : inline void pup( PUP::er& p, Term& t ) { t.pup(p); }
93 : :
94 : : //! \brief Products are arbitrary number of Terms to be multiplied and ensemble
95 : : //! averaged.
96 : : //! \details An example is the scalar flux in x direction which needs two terms
97 : : //! for ensemble averaging: (Y-\<Y\>) and (U-\<U\>), then the central moment is
98 : : //! \<yu\> = <(Y-\<Y\>)(U-\<U\>)>, another example is the third mixed central
99 : : //! moment of three scalars which needs three terms for ensemble averaging:
100 : : //! (Y1-\<Y1\>), (Y2-\<Y2\>), and (Y3-\<Y3\>), then the central moment is
101 : : //! \<y1y2y3\> = \<(Y1-\<Y1\>)(Y2-\<Y2\>)(Y3-\<Y3\>)\>.
102 : : using Product = std::vector< Term >;
103 : :
104 : : //! \brief Case-insensitive character comparison functor
105 : : struct CaseInsensitiveCharLess {
106 : : //! Function call operator
107 : : //! \param[in] lhs Left character of the comparitor operand
108 : : //! \param[in] rhs Right character of the comparitor operand
109 : : //! \return Boolean indicating the result of the comparison
110 : : bool operator() ( char lhs, char rhs ) const {
111 : : return std::tolower( lhs ) < std::tolower( rhs );
112 : : }
113 : : };
114 : :
115 : : //! \brief Find out if a vector of Terms only contains ordinary moment terms
116 : : //! \details If and only if all terms are ordinary, the vector of Terms is
117 : : //! ordinary.
118 : : //! \param[in] vec Vector of Terms to check
119 : : //! \return Boolean indicating if all terms are ordinary
120 : : static inline bool
121 : : ordinary( const std::vector< ctr::Term >& vec ) {
122 : : if (std::any_of( vec.cbegin(), vec.cend(),
123 : 0 : []( const ctr::Term& t ){ return t.moment == ctr::Moment::CENTRAL; } ))
124 : : return false;
125 : : else
126 : : return true;
127 : : }
128 : :
129 : : //! \brief Find out if a vector of Terms contains any central moment terms
130 : : //! \details If any term is central, the vector of Terms is central.
131 : : //! \param[in] vec Vector of Terms to check
132 : : //! \return Boolean indicating of any term is central
133 : : static inline bool
134 : : central( const std::vector< ctr::Term >& vec )
135 : : { return !ordinary( vec ); }
136 : :
137 : : //! \brief Probability density function (vector of sample space variables)
138 : : using Probability = std::vector< Term >;
139 : :
140 : : //! \brief PDF information bundle
141 : : //! \note If the user did not specify extents for a PDF, the corresponding
142 : : //! extents vector still exists but it is empty.
143 : : struct PDFInfo {
144 : : const std::string& name; //!< PDF identifier, i.e., name
145 : : const std::vector< tk::real >& exts; //!< Sample space extents
146 : : std::vector< std::string > vars; //!< List of sample space variables
147 : : std::uint64_t it; //!< Iteration count
148 : : tk::real time; //!< Time stamp
149 : : };
150 : :
151 : : //! \brief Find PDF information, see tk::ctr::PDFInfo
152 : : //! \note Size of binsizes, names, pdfs, and exts must all be equal
153 : : //! \note idx must be less than the length of binsizes, names, and pdfs
154 : : //! \param[in] binsizes Vector of sample space bin sizes for multiple PDFs
155 : : //! \param[in] names Vector of PDF names
156 : : //! \param[in] exts Vector of sample space extents. Note: if the user did not
157 : : //! specify extents for a PDF, the corresponding extents vector still exists
158 : : //! but it is empty.
159 : : //! \param[in] pdfs Vector of PDFs
160 : : //! \param[in] m ORDINARY or CENTRAL PDF we are looking for
161 : : //! \param[in] idx Index of the PDF within the list of matching (based on D and
162 : : //! m) PDFs requested
163 : : //! \param[in] it Iteration count
164 : : //! \param[in] time Physical time
165 : : //! \return The PDF metadata requested
166 : : //! \details Find PDF information given the sample space dimension (template
167 : : //! argument D), ordinary or central PDF (m), and the index within the list of
168 : : //! matching (based on D and m) PDFs requested (idx). This function must find
169 : : //! the PDF, if it does not, it throws an exception.
170 : : //! \see walker::Distributor
171 : : template< std::size_t D >
172 : : PDFInfo pdfInfo( const std::vector< std::vector< tk::real > >& binsizes,
173 : : const std::vector< std::string >& names,
174 : : const std::vector< std::vector< tk::real > >& exts,
175 : : const std::vector< Probability >& pdfs,
176 : : tk::ctr::Moment m,
177 : : std::size_t idx,
178 : : std::uint64_t it,
179 : : tk::real time )
180 : : {
181 : : Assert( binsizes.size() == names.size(),
182 : : "Length of binsizes vector and that of PDF names must equal" );
183 : : Assert( binsizes.size() == pdfs.size(),
184 : : "Length of binsizes vector and that of PDFs must equal" );
185 : : Assert( binsizes.size() == exts.size(),
186 : : "Length of binsizes vector and that of PDF extents must equal" );
187 : : Assert( binsizes.size() > idx, "Indexing out of bounds" );
188 : :
189 : : std::size_t i = 0; // will count all PDFs queried
190 : : std::size_t n = 0; // will count PDFs with sample space dimensions and type
191 : : // (orindary or central) we are looking for
192 : : for (const auto& bs : binsizes) {
193 : : if ( bs.size() == D &&
194 : : ((m == Moment::ORDINARY && ordinary(pdfs[i])) ||
195 : : (m == Moment::CENTRAL && central(pdfs[i]))) ) ++n;
196 : : if (n == idx+1) {
197 : : std::vector< std::string > vars;
198 : : for (const auto& term : pdfs[i])
199 : : // cppcheck-suppress useStlAlgorithm
200 : : vars.push_back( term.var + std::to_string(term.field+1) );
201 : : return { names[i], exts[i], std::move(vars), it, time };
202 : : }
203 : : ++i;
204 : : }
205 : : Throw( "Cannot find PDF." );
206 : : }
207 : :
208 : : //! Construct mean
209 : : //! \param[in] var Variable
210 : : //! \param[in] c Component number
211 : : //! \return Constructed vector< Term > identifying the first ordinary moment
212 : : //! (mean) of field (component) c of variable var
213 : : static inline Product
214 : 4 : mean( char var, uint64_t c ) {
215 : 4 : tk::ctr::Term m( static_cast<char>(std::toupper(var)), c, Moment::ORDINARY );
216 [ + - ]: 4 : return tk::ctr::Product( { m } );
217 : : }
218 : :
219 : : //! Construct variance
220 : : //! \param[in] var Variable
221 : : //! \param[in] c Component number
222 : : //! \return Constructed vector< Term > identifying the second central moment
223 : : //! (variance) of field (component) c of variable var
224 : : static inline Product
225 : 8 : variance( char var, uint64_t c ) {
226 : 8 : tk::ctr::Term f( static_cast<char>(std::tolower(var)), c, Moment::CENTRAL );
227 [ + - ]: 8 : return tk::ctr::Product( { f, f } );
228 : : }
229 : :
230 : : } // ctr::
231 : : } // tk::
232 : :
233 : : #endif // StatControl_h
|
