Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/FluxCorrector.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 FluxCorrector performs limiting for transport equations
9 : : \details FluxCorrector performs limiting for transport equations. Each
10 : : FluxCorrector object performs the limiting procedure, according to a
11 : : flux-corrected transport algorithm, on a chunk of the full load (part of the
12 : : mesh).
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef FluxCorrector_h
16 : : #define FluxCorrector_h
17 : :
18 : : #include <utility>
19 : : #include <unordered_map>
20 : :
21 : : #include "NoWarning/pup.hpp"
22 : :
23 : : #include "Keywords.hpp"
24 : : #include "Fields.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : :
27 : : namespace inciter {
28 : :
29 : : extern ctr::InputDeck g_inputdeck;
30 : :
31 : : //! FluxCorrector is used to perform flux-corrected transport
32 : : //! \see Löhner, R., Morgan, K., Peraire, J. and Vahdati, M. (1987), Finite
33 : : //! element flux-corrected transport (FEM–FCT) for the Euler and Navier–Stokes
34 : : //! equations. Int. J. Numer. Meth. Fluids, 7: 1093–1109.
35 : : //! doi:10.1002/fld.1650071007
36 : : class FluxCorrector {
37 : :
38 : : private:
39 : : using ncomp_t = kw::ncomp::info::expect::type;
40 : :
41 : : public:
42 : : //! Constructor
43 : : //! \param[in] is Size of the mesh element connectivity vector (inpoel size)
44 : 2576 : explicit FluxCorrector( std::size_t is = 0 ) :
45 : 2576 : m_aec( is, g_inputdeck.get< tag::component >().nprop() ),
46 : : m_sys( findsys< tag::compflow >() ),
47 [ + - ][ + - ]: 2576 : m_vel( findvel< tag::compflow >() ) {}
48 : :
49 : : //! Collect scalar comonent indices for equation systems
50 : : //! \tparam Eq Equation types to consider as equation systems
51 : : //! \return List of component indices to treat as a system
52 : : template< class... Eq >
53 : : std::vector< std::vector< ncomp_t > >
54 : 2576 : findsys() {
55 : 2576 : std::vector< std::vector< ncomp_t > > sys;
56 : 2576 : ( ... , [&](){
57 : : // Access system-FCT variable indices for all systems of type Eq
58 : 2576 : const auto& sv = g_inputdeck.get< tag::param, Eq, tag::sysfctvar >();
59 : : // Access number of scalar components in all systems of type Eq
60 : 2576 : const auto& ncompv = g_inputdeck.get< tag::component >().get< Eq >();
61 : : // Assign variable indices for system FCT for each Eq system
62 [ + + ]: 3401 : for (std::size_t e=0; e<ncompv.size(); ++e) {
63 [ + + ]: 825 : if (g_inputdeck.get< tag::param, Eq, tag::sysfct >().at(e)) {
64 : 96 : auto offset = g_inputdeck.get< tag::component >().offset< Eq >( e );
65 [ + - ]: 96 : sys.push_back( std::vector< ncomp_t >() );
66 [ + - ][ + + ]: 576 : for (auto c : sv.at(e)) {
67 [ + - ]: 480 : sys.back().push_back( offset + c );
68 : : }
69 : : }
70 [ + - ]: 2576 : } }() );
71 [ + + ]: 2672 : for ([[maybe_unused]] const auto& s : sys) {
72 [ + - ][ - + ]: 576 : Assert( std::all_of( begin(s), end(s), [&]( std::size_t i ){
[ - - ][ - - ]
[ - - ]
73 : : return i < g_inputdeck.get< tag::component >().nprop(); } ),
74 : : "Eq system index larger than total number of components" );
75 : : }
76 : 2576 : return sys;
77 : : }
78 : :
79 : : //! Find components of a velocity for equation systems
80 : : //! \tparam Eq Equation types to consider as equation systems
81 : : //! \return List of 3 component indices to treat as a velocity
82 : : //! \warning Currently, this is only a punt for single-material flow: we
83 : : //! simply take the components 1,2,3 as the velocity for each system of
84 : : //! type Eq
85 : : template< class... Eq >
86 : : std::vector< std::array< ncomp_t, 3 > >
87 : 2576 : findvel() {
88 : 2576 : std::vector< std::array< ncomp_t, 3 > > vel;
89 : 2576 : ( ... , [&](){
90 : : // Access number of scalar components in all systems of type Eq
91 : 2576 : const auto& ncompv = g_inputdeck.get< tag::component >().get< Eq >();
92 : : // Assign variable indices for system FCT for each Eq system
93 [ + + ]: 3401 : for (std::size_t e=0; e<ncompv.size(); ++e) {
94 : 825 : auto offset = g_inputdeck.get< tag::component >().offset< Eq >( e );
95 [ + - ]: 825 : vel.push_back( { offset+1, offset+2, offset+3 } );
96 [ + - ]: 2576 : } }() );
97 [ + + ]: 3401 : for ([[maybe_unused]] const auto& v : vel) {
98 [ + - ][ - + ]: 3300 : Assert( std::all_of( begin(v), end(v), [&]( std::size_t i ){
[ - - ][ - - ]
[ - - ]
99 : : return i < g_inputdeck.get< tag::component >().nprop(); } ),
100 : : "Eq system index larger than total number of components" );
101 : : }
102 : 2576 : return vel;
103 : : }
104 : :
105 : : //! Resize state (e.g., after mesh refinement)
106 : 76 : void resize( std::size_t is ) { m_aec.resize( is ); }
107 : :
108 : : //! Compute antidiffusive element contributions (AEC)
109 : : void aec(
110 : : const std::array< std::vector< tk::real >, 3 >& coord,
111 : : const std::vector< std::size_t >& inpoel,
112 : : const std::vector< tk::real >& vol,
113 : : const std::unordered_map< std::size_t,
114 : : std::vector< std::pair< bool, tk::real > > >& bcdir,
115 : : const std::unordered_map< int,
116 : : std::unordered_set< std::size_t > >& symbcnodemap,
117 : : const std::unordered_map< int,
118 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
119 : : const tk::Fields& Un,
120 : : tk::Fields& P );
121 : :
122 : : //! Verify the assembled antidiffusive element contributions
123 : : bool verify( std::size_t nchare,
124 : : const std::vector< std::size_t >& inpoel,
125 : : const tk::Fields& dUh,
126 : : const tk::Fields& dUl ) const;
127 : :
128 : : //! Compute mass diffusion contribution to the rhs of the low order system
129 : : tk::Fields diff( const std::array< std::vector< tk::real >, 3 >& coord,
130 : : const std::vector< std::size_t >& inpoel,
131 : : const tk::Fields& Un ) const;
132 : :
133 : : //! \brief Compute the maximum and minimum unknowns of all elements
134 : : //! surrounding nodes
135 : : void alw( const std::vector< std::size_t >& inpoel,
136 : : const tk::Fields& Un,
137 : : const tk::Fields& Ul,
138 : : tk::Fields& Q ) const;
139 : :
140 : : //! Compute limited antiffusive element contributions and apply to mesh nodes
141 : : void lim( const std::vector< std::size_t >& inpoel,
142 : : const std::unordered_map< std::size_t,
143 : : std::vector< std::pair< bool, tk::real > > >& bcdir,
144 : : const tk::Fields& P,
145 : : const tk::Fields& Ul,
146 : : tk::Fields& Q,
147 : : tk::Fields& A ) const;
148 : :
149 : : // Collect mesh output fields from FCT
150 : : std::tuple< std::vector< std::string >,
151 : : std::vector< std::vector< tk::real > > >
152 : : fields( const std::vector< std::size_t >& inpoel ) const;
153 : :
154 : : /** @name Charm++ pack/unpack serializer member functions */
155 : : ///@{
156 : : //! \brief Pack/Unpack serialize member function
157 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
158 : 5766 : void pup( PUP::er& p ) {
159 : 5766 : p | m_aec;
160 : 5766 : p | m_sys;
161 : 5766 : p | m_vel;
162 : 5766 : }
163 : : //! \brief Pack/Unpack serialize operator|
164 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
165 : : //! \param[in,out] i FluxCorrector object reference
166 : 5766 : friend void operator|( PUP::er& p, FluxCorrector& i ) { i.pup(p); }
167 : : //@}
168 : :
169 : : private:
170 : : //! Antidiffusive element contributions for all scalar components
171 : : tk::Fields m_aec;
172 : : //! Component indices to treat as a system for multiple systems
173 : : std::vector< std::vector< ncomp_t > > m_sys;
174 : : //! Component indices to treat as a velocity vector for multiple systems
175 : : std::vector< std::array< ncomp_t, 3 > > m_vel;
176 : : };
177 : :
178 : : } // inciter::
179 : :
180 : : #endif // FluxCorrector_h
|
