Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/FluxCorrector.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 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 : :
16 : : #include <limits>
17 : : #include <sstream>
18 : : #include <algorithm>
19 : : #include <unordered_map>
20 : :
21 : : #include "Macro.hpp"
22 : : #include "Vector.hpp"
23 : : #include "Around.hpp"
24 : : #include "DerivedData.hpp"
25 : : #include "FluxCorrector.hpp"
26 : : #include "Inciter/InputDeck/InputDeck.hpp"
27 : :
28 : : using inciter::FluxCorrector;
29 : :
30 : : void
31 : 18283 : FluxCorrector::aec(
32 : : const std::array< std::vector< tk::real >, 3 >& coord,
33 : : const std::vector< std::size_t >& inpoel,
34 : : const std::vector< tk::real >& vol,
35 : : const std::unordered_map< std::size_t,
36 : : std::vector< std::pair< bool, tk::real > > >& bcdir,
37 : : const std::unordered_map< int,
38 : : std::unordered_set< std::size_t > >& symbcnodemap,
39 : : const std::unordered_map< int,
40 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
41 : : const tk::Fields& Un,
42 : : tk::Fields& P )
43 : : // *****************************************************************************
44 : : // Compute antidiffusive element contributions (AEC)
45 : : //! \param[in] coord Mesh node coordinates
46 : : //! \param[in] inpoel Mesh element connectivity
47 : : //! \param[in] vol Volume associated to mesh nodes
48 : : //! \param[in] bcdir Vector of pairs of bool and boundary condition value
49 : : //! associated to local mesh node IDs at which to set Dirichlet boundary
50 : : //! conditions.
51 : : //! \param[in] symbcnodemap Unique set of node ids at which to set symmetry BCs
52 : : //! associated to side set ids
53 : : //! \param[in] bnorm Face normals in boundary points: key global node id,
54 : : //! value: unit normal, outer key: side set id
55 : : //! \param[in] Un Solution at the previous time step
56 : : //! \param[in,out] P The sums of positive (negative) AECs to nodes
57 : : //! \details The antidiffusive element contributions (AEC) are defined as the
58 : : //! difference between the high and low order solution, where the high order
59 : : //! solution is obtained from consistent mass Taylor-Galerkin discretization
60 : : //! and the low order solution is lumped mass Taylor-Galerkin + diffusion.
61 : : //! Note that AEC is not directly computed as dUh - dUl (although that could
62 : : //! also be done), but as AEC = M_L^{-1} (M_Le - M_ce) (ctau * Un + dUh),
63 : : //! where
64 : : //! * M_ce is the element's consistent mass matrix,
65 : : //! * M_Le is the element's lumped mass matrix,
66 : : //! * ctau is the mass diffusion coefficient on the rhs of the low order
67 : : //! solution, see also FluxCorrector::diff(),
68 : : //! * Un is the solution at the previous time step
69 : : //! * dUh is the increment of the high order solution, and
70 : : //! * M_L^{-1} is the inverse of the assembled lumped mass matrix, i.e., the
71 : : //! volume associated to a mesh node by summing the quarter of the element
72 : : //! volumes surrounding the node. Note that this is the correct node volume
73 : : //! taking into account that some nodes are on chare boundaries.
74 : : //! \note Since we use the lumped-mass for the high-order solution, dUh
75 : : //! does not contribute to AEC, as computed above.
76 : : //! \see Löhner, R., Morgan, K., Peraire, J. and Vahdati, M. (1987), Finite
77 : : //! element flux-corrected transport (FEM–FCT) for the Euler and Navier–Stokes
78 : : //! equations. Int. J. Numer. Meth. Fluids, 7: 1093–1109.
79 : : //! doi:10.1002/fld.1650071007
80 : : // *****************************************************************************
81 : : {
82 : 18283 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
83 : 18283 : auto ctau = g_inputdeck.get< tag::discr, tag::ctau >();
84 : :
85 : : Assert( vol.size() == coord[0].size(), "Nodal volume vector size mismatch" );
86 : : Assert( m_aec.nunk() == inpoel.size() && m_aec.nprop() == ncomp,
87 : : "AEC and mesh connectivity size mismatch" );
88 : : Assert( Un.nunk() == P.nunk() && Un.nprop() == P.nprop()/2, "Size mismatch" );
89 : :
90 : : const auto& x = coord[0];
91 : : const auto& y = coord[1];
92 : : const auto& z = coord[2];
93 : :
94 : : m_aec.fill( 0.0 );
95 : :
96 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
97 : 9254852 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
98 : 9254852 : inpoel[e*4+2], inpoel[e*4+3] }};
99 : :
100 : : // compute element Jacobi determinant
101 : : const std::array< tk::real, 3 >
102 : 9254852 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
103 : 9254852 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
104 : 9254852 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
105 : : const auto J = tk::triple( ba, ca, da );
106 : : Assert( J > 0, "Element Jacobian non-positive" );
107 : :
108 : : // lumped - consistent mass
109 : : std::array< std::array< tk::real, 4 >, 4 > m; // nnode*nnode [4][4]
110 : 9254852 : m[0][0] = m[1][1] = m[2][2] = m[3][3] = 3.0*J/120.0;// diagonal
111 : 9254852 : m[0][1] = m[0][2] = m[0][3] = // off-diagonal
112 : 9254852 : m[1][0] = m[1][2] = m[1][3] =
113 : 9254852 : m[2][0] = m[2][1] = m[2][3] =
114 : 9254852 : m[3][0] = m[3][1] = m[3][2] = -J/120.0;
115 : :
116 : : // access solution at element nodes at time n
117 : 9254852 : std::vector< std::array< tk::real, 4 > > un( ncomp );
118 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) un[c] = Un.extract( c, 0, N );
119 : :
120 : : // Compute antidiffusive element contributions (AEC). The high order system
121 : : // is M_c * dUh = r, where M_c is the consistent mass matrix and r is the
122 : : // high order right hand side. The low order system is constructed from the
123 : : // high order one by lumping the consistent mass matrix and adding mass
124 : : // diffusion: M_L * dUl = r + c_tau * (M_c - M_L) Un, where M_L is the
125 : : // lumped mass matrix, c_tau is the mass diffusion coefficient (c_tau = 1.0
126 : : // guarantees a monotonic solution). See also the details in the function
127 : : // header for the notation. Based on the above, the AEC, in general, is
128 : : // computed as AEC = M_L^{-1} (M_Le - M_ce) (ctau * Un + dUh), which can be
129 : : // obtained by subtracting the low order system from the high order system.
130 : : // Note that the solution update is U^{n+1} = Un + dUl + lim(dUh - dUl),
131 : : // where the last term is the limited AEC. (Think of 'lim' as the limit
132 : : // coefficient between 0 and 1.)
133 [ + + ]: 46274260 : for (std::size_t j=0; j<4; ++j)
134 [ + + ]: 92573776 : for (ncomp_t c=0; c<ncomp; ++c)
135 [ + + ]: 277771840 : for (std::size_t k=0; k<4; ++k)
136 : 222217472 : m_aec(e*4+j,c,0) += m[j][k] * ctau*un[c][k] / vol[N[j]];
137 : : }
138 : :
139 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
140 : : const std::array< std::size_t, 4 >
141 : 9254852 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
142 [ + + ]: 46274260 : for (std::size_t j=0; j<4; ++j) {
143 : : // Dirichlet BCs: At nodes where Dirichlet boundary conditions (BC) are
144 : : // set, we set the AEC to zero. This is because if the (same) BCs are
145 : : // correctly set for both the low and the high order solution, there
146 : : // should be no difference between the low and high order increments,
147 : : // thus AEC = dUh - dUl = 0.
148 : : auto b = bcdir.find(N[j]);
149 [ + + ]: 37019408 : if (b != end(bcdir)) {
150 [ + + ]: 18417696 : for (ncomp_t c=0; c<ncomp; ++c) {
151 [ + - ]: 14546180 : if (b->second[c].first) {
152 : 14546180 : m_aec(e*4+j,c,0) = 0.0;
153 : : }
154 : : }
155 : : }
156 : : // Symmetry BCs
157 [ + + ]: 37787488 : for (const auto& [s,nodes] : symbcnodemap) {
158 : : auto i = nodes.find(N[j]);
159 [ + + ]: 768080 : if (i != end(nodes)) {
160 : : auto l = bnorm.find(s);
161 [ + - ]: 184280 : if (l != end(bnorm)) {
162 : : auto k = l->second.find(N[j]);
163 [ + + ]: 368560 : for (const auto& vel : m_vel) {
164 : : std::array< tk::real, 3 >
165 : 184280 : v{ m_aec(e*4+j,vel[0],0),
166 : 184280 : m_aec(e*4+j,vel[1],0),
167 : 184280 : m_aec(e*4+j,vel[2],0) },
168 : 184280 : n{ k->second[0], k->second[1], k->second[2] };
169 : : auto vn = tk::dot( v, n );
170 : 184280 : m_aec(e*4+j,vel[0],0) -= vn * n[0];
171 : 184280 : m_aec(e*4+j,vel[1],0) -= vn * n[1];
172 : 184280 : m_aec(e*4+j,vel[2],0) -= vn * n[2];
173 : : }
174 : : }
175 : : }
176 : : }
177 : : }
178 : : }
179 : :
180 : : // sum all positive (negative) antidiffusive element contributions to nodes
181 : : // (Lohner: P^{+,-}_i)
182 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
183 : 9254852 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
184 : 9254852 : inpoel[e*4+2], inpoel[e*4+3] }};
185 [ + + ]: 46274260 : for (std::size_t j=0; j<4; ++j) {
186 [ + + ]: 92573776 : for (ncomp_t c=0; c<ncomp; ++c) {
187 [ + + ][ + + ]: 69822841 : P(N[j],c*2+0,0) += std::max( 0.0, m_aec(e*4+j,c,0) );
188 [ + + ]: 76153866 : P(N[j],c*2+1,0) += std::min( 0.0, m_aec(e*4+j,c,0) );
189 : : }
190 : : }
191 : : }
192 : 18283 : }
193 : :
194 : : bool
195 : 18283 : FluxCorrector::verify( std::size_t nchare,
196 : : const std::vector< std::size_t >& inpoel,
197 : : const tk::Fields& dUh,
198 : : const tk::Fields& dUl ) const
199 : : // *****************************************************************************
200 : : // Verify the assembled antidiffusive element contributions (AEC)
201 : : //! \param[in] nchare Total number of host chares
202 : : //! \param[in] inpoel Mesh element connectivity
203 : : //! \param[in] dUh Increment of the high order solution
204 : : //! \param[in] dUl Increment of the low order solution
205 : : //! \return True if verification has been done
206 : : //! \details This verification only makes sense if no communication is to be
207 : : //! done, i.e., if there is a single host chare, because the AEC assembled to
208 : : //! mesh nodes only contains partial contributions on chare boundaries at this
209 : : //! point. Verification in parallel would incure communication of the
210 : : //! unlimited AEC, which in general is not necessary, so we will not do that
211 : : //! for the sake of verification.
212 : : //! \note Client code should ensure that this function is optimized away in
213 : : //! RELEASE mode.
214 : : // *****************************************************************************
215 : : {
216 : : Assert( dUl.nunk() == dUh.nunk() && dUl.nprop() == dUh.nprop(),
217 : : "Unknown array size mismatch" );
218 : :
219 [ + + ]: 18283 : if (nchare == 1) {
220 : 499 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
221 : 499 : tk::Fields U( dUh.nunk(), dUh.nprop() );
222 : : U.fill( 0.0 );
223 : :
224 [ + + ]: 1989253 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
225 [ + - ]: 1988754 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
226 : 1988754 : inpoel[e*4+2], inpoel[e*4+3] }};
227 : : // access pointer to solution at element nodes
228 [ + - ][ - - ]: 1988754 : std::vector< const tk::real* > u( ncomp );
229 [ + + ]: 5251688 : for (ncomp_t c=0; c<ncomp; ++c) u[c] = U.cptr( c, 0 );
230 : : // scatter-add antidiffusive element contributions to nodes
231 [ + + ]: 5251688 : for (ncomp_t c=0; c<ncomp; ++c)
232 [ + + ]: 16314670 : for (std::size_t j=0; j<4; ++j)
233 : 13051736 : U.var(u[c],N[j]) += m_aec(e*4+j,c,0);
234 : : }
235 : :
236 : : // Compute maximum difference between the assembled AEC and dUh-dUl
237 [ + - ][ + - ]: 499 : auto d = tk::maxdiff( U, dUh-dUl );
238 : :
239 : : // Tolerance: 10 x the linear solver tolerance for the high order solution.
240 [ - + ]: 499 : if (d.second > 1.0e-7) {
241 : : const auto& duh = dUh.vec();
242 : : const auto& dul = dUl.vec();
243 : : const auto& u = U.vec();
244 [ - - ]: 0 : std::stringstream ss;
245 : : ss << "maximum difference at mesh node " << d.first << ": " << d.second
246 [ - - ][ - - ]: 0 : << ", dUh:" << duh[d.first] << ", dUl:" << dul[d.first]
247 [ - - ]: 0 : << ", dUh-dUl:" << duh[d.first] - dul[d.first]
248 [ - - ]: 0 : << ", AEC:" << u[d.first];
249 [ - - ][ - - ]: 0 : Throw( "Assembled AEC does not equal dUh-dUl: " + ss.str() );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
250 : : }
251 : :
252 : : return true;
253 : : }
254 : :
255 : : return false;
256 : : }
257 : :
258 : : tk::Fields
259 : 18283 : FluxCorrector::diff( const std::array< std::vector< tk::real >, 3 >& coord,
260 : : const std::vector< std::size_t >& inpoel,
261 : : const tk::Fields& Un ) const
262 : : // *****************************************************************************
263 : : // Compute mass diffusion contribution to the RHS of the low order system
264 : : //! \param[in] coord Mesh node coordinates
265 : : //! \param[in] inpoel Mesh element connectivity
266 : : //! \param[in] Un Solution at the previous time step
267 : : //! \return Mass diffusion contribution to the RHS of the low order system
268 : : // *****************************************************************************
269 : : {
270 : 18283 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
271 : 18283 : auto ctau = g_inputdeck.get< tag::discr, tag::ctau >();
272 : :
273 : : // access node coordinates
274 : : const auto& x = coord[0];
275 : : const auto& y = coord[1];
276 : : const auto& z = coord[2];
277 : :
278 : 18283 : tk::Fields D( Un.nunk(), Un.nprop() );
279 : : D.fill( 0.0 );
280 : :
281 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
282 : : // access node IDs
283 : : const std::array< std::size_t, 4 >
284 [ + - ]: 9254852 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
285 : : // compute element Jacobi determinant
286 : : const std::array< tk::real, 3 >
287 : 9254852 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
288 : 9254852 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
289 [ + - ]: 9254852 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
290 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
291 : : Assert( J > 0, "Element Jacobian non-positive" );
292 : :
293 : : // lumped - consistent mass
294 : : std::array< std::array< tk::real, 4 >, 4 > m; // nnode*nnode [4][4]
295 : 9254852 : m[0][0] = m[1][1] = m[2][2] = m[3][3] = 3.0*J/120.0;// diagonal
296 : 9254852 : m[0][1] = m[0][2] = m[0][3] = // off-diagonal
297 : 9254852 : m[1][0] = m[1][2] = m[1][3] =
298 : 9254852 : m[2][0] = m[2][1] = m[2][3] =
299 : 9254852 : m[3][0] = m[3][1] = m[3][2] = -J/120.0;
300 : :
301 : : // access solution at element nodes at time n
302 [ + - ][ - - ]: 9254852 : std::vector< std::array< tk::real, 4 > > un( ncomp );
303 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) un[c] = Un.extract( c, 0, N );
304 : : // access pointer to mass diffusion right hand side at element nodes
305 [ + - ][ - - ]: 9254852 : std::vector< const tk::real* > d( ncomp );
306 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) d[c] = D.cptr( c, 0 );
307 : :
308 : : // scatter-add mass diffusion element contributions to rhs nodes
309 [ + + ]: 46274260 : for (std::size_t a=0; a<4; ++a) {
310 [ + + ]: 92573776 : for (ncomp_t c=0; c<ncomp; ++c)
311 [ + + ]: 277771840 : for (std::size_t b=0; b<4; ++b)
312 : 222217472 : D.var(d[c],N[a]) -= ctau * m[a][b] * un[c][b];
313 : : }
314 : : }
315 : :
316 : 18283 : return D;
317 : : }
318 : :
319 : : void
320 : 18283 : FluxCorrector::alw( const std::vector< std::size_t >& inpoel,
321 : : const tk::Fields& Un,
322 : : const tk::Fields& Ul,
323 : : tk::Fields& Q ) const
324 : : // *****************************************************************************
325 : : // Compute the maximum and minimum unknowns of elements surrounding nodes
326 : : //! \param[in] inpoel Mesh element connectivity
327 : : //! \param[in] Un Solution at the previous time step
328 : : //! \param[in] Ul Low order solution
329 : : //! \param[in,out] Q Maximum and mimimum unknowns of elements surrounding nodes
330 : : // *****************************************************************************
331 : : {
332 : : Assert( Q.nunk() == Un.nunk() && Q.nprop() == Un.nprop()*2, "Max and min "
333 : : "unknowns of elements surrounding nodes array size mismatch" );
334 : :
335 : 18283 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
336 : 18283 : auto clip = g_inputdeck.get< tag::discr, tag::fctclip >();
337 : :
338 : : // compute maximum and minimum nodal values of all elements (Lohner: u^*_el)
339 : 18283 : tk::Fields S( inpoel.size()/4, ncomp*2 );
340 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
341 : 9254852 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
342 : 9254852 : inpoel[e*4+2], inpoel[e*4+3] }};
343 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) {
344 : 13888592 : S(e,c*2+0,0) = -std::numeric_limits< tk::real >::max();
345 : 13888592 : S(e,c*2+1,0) = std::numeric_limits< tk::real >::max();
346 [ + + ]: 69442960 : for (std::size_t j=0; j<4; ++j) {
347 : : // compute maximum and minimum nodal values of Ul and Un (Lohner: u^*_i)
348 [ - + ][ + + ]: 64957178 : auto jmax = clip ? Ul(N[j],c,0) : std::max(Ul(N[j],c,0), Un(N[j],c,0));
349 [ - + ][ + + ]: 85038849 : auto jmin = clip ? Ul(N[j],c,0) : std::min(Ul(N[j],c,0), Un(N[j],c,0));
350 [ + + ]: 55554368 : if (jmax > S(e,c*2+0,0)) S(e,c*2+0,0) = jmax;
351 [ + + ]: 55554368 : if (jmin < S(e,c*2+1,0)) S(e,c*2+1,0) = jmin;
352 : : }
353 : : }
354 : : }
355 : :
356 : : // compute maximum and mimimum unknowns of all elements surrounding each node
357 : : // (Lohner: u^{max,min}_i)
358 [ + - ]: 36566 : const auto esup = tk::genEsup( inpoel, 4 );
359 [ + + ]: 2755420 : for (std::size_t p=0; p<Un.nunk(); ++p) {
360 [ + + ]: 39756545 : for (auto e : tk::Around(esup,p)) {
361 [ + + ]: 92573776 : for (ncomp_t c=0; c<ncomp; ++c) {
362 [ + + ]: 55554368 : if (S(e,c*2+0,0) > Q(p,c*2+0,0)) Q(p,c*2+0,0) = S(e,c*2+0,0);
363 [ + + ]: 55554368 : if (S(e,c*2+1,0) < Q(p,c*2+1,0)) Q(p,c*2+1,0) = S(e,c*2+1,0);
364 : : }
365 : : }
366 : : }
367 : 18283 : }
368 : :
369 : : void
370 : 18283 : FluxCorrector::lim( const std::vector< std::size_t >& inpoel,
371 : : const std::unordered_map< std::size_t,
372 : : std::vector< std::pair< bool, tk::real > > >& bcdir,
373 : : const tk::Fields& P,
374 : : const tk::Fields& Ul,
375 : : tk::Fields& Q,
376 : : tk::Fields& A ) const
377 : : // *****************************************************************************
378 : : // Compute limited antiffusive element contributions and apply to mesh nodes
379 : : //! \param[in] inpoel Mesh element connectivity
380 : : //! \param[in] bcdir Vector of pairs of bool and boundary condition value
381 : : //! associated to mesh node IDs at which to set Dirichlet boundary conditions.
382 : : //! \param[in] P The sums of all positive (negative) AECs to nodes
383 : : //! \param[in] Ul Low order solution
384 : : //! \param[in,out] Q The maximum and mimimum unknowns of elements surrounding
385 : : //! each node
386 : : //! \param[in,out] A Limited antidiffusive element contributions scatter-added
387 : : //! to nodes
388 : : //! \note Q is also overwritten to avoid using temporary memory
389 : : // *****************************************************************************
390 : : {
391 : : Assert( P.nunk() == Q.nunk() && P.nprop() == Q.nprop(), "Size mismatch" );
392 : : Assert( P.nunk() == Ul.nunk() && P.nprop() == Ul.nprop()*2, "Size mismatch" );
393 : : Assert( A.nunk() == Ul.nunk() && A.nprop() == Ul.nprop(), "Size mismatch" );
394 : :
395 : 18283 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
396 : :
397 : : // compute the maximum and minimum increments and decrements nodal solution
398 : : // values are allowed to achieve (Lohner: Q^{+,-}_i)
399 [ + + ]: 2755420 : for (std::size_t p=0; p<Ul.nunk(); ++p)
400 [ + + ]: 7596754 : for (ncomp_t c=0; c<ncomp; ++c) {
401 : 4859617 : Q(p,c*2+0,0) -= Ul(p,c,0);
402 : 4859617 : Q(p,c*2+1,0) -= Ul(p,c,0);
403 : : }
404 : :
405 : 18283 : auto eps = g_inputdeck.get< tag::discr, tag::fcteps >();
406 : :
407 : : // compute the ratios of positive and negative element contributions that
408 : : // ensure monotonicity (Lohner: R^{+,-})
409 [ + + ]: 2755420 : for (std::size_t p=0; p<P.nunk(); ++p) {
410 [ + + ]: 7596754 : for (ncomp_t c=0; c<ncomp; ++c) {
411 : :
412 [ + + ]: 4859617 : if (P(p,c*2+0,0) < eps)
413 : 3077098 : Q(p,c*2+0,0) = 1.0;
414 : : else
415 [ + + ]: 1848880 : Q(p,c*2+0,0) = std::min(1.0,Q(p,c*2+0,0)/P(p,c*2+0,0));
416 : :
417 [ + + ]: 4859617 : if (P(p,c*2+1,0) > -eps)
418 : 2665865 : Q(p,c*2+1,0) = 1.0;
419 : : else
420 [ + + ]: 2652474 : Q(p,c*2+1,0) = std::min(1.0,Q(p,c*2+1,0)/P(p,c*2+1,0));
421 : :
422 : : }
423 : : }
424 : :
425 : : // calculate limit coefficient for all elements (Lohner: C_el)
426 : 18283 : tk::Fields C( inpoel.size()/4, ncomp );
427 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
428 : 9254852 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
429 : 9254852 : inpoel[e*4+2], inpoel[e*4+3] }};
430 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) {
431 : : std::array< tk::real, 4 > R;
432 [ + + ]: 69442960 : for (std::size_t j=0; j<4; ++j) {
433 : :
434 [ + + ]: 55554368 : if (std::abs(m_aec(e*4+j,c,0)) < eps)
435 : 23434102 : R[j] = 1.0;
436 [ + + ]: 32120266 : else if (m_aec(e*4+j,c,0) > 0.0)
437 : 13060745 : R[j] = Q(N[j],c*2+0,0);
438 : : else
439 : 19059521 : R[j] = Q(N[j],c*2+1,0);
440 : :
441 : : }
442 [ - + ]: 13888592 : C(e,c,0) = *std::min_element( begin(R), end(R) );
443 : : // if all vertices happened to be on a Dirichlet boundary, ignore limiting
444 [ - + ]: 13888592 : if (C(e,c,0) > 1.0) C(e,c,0) = 1.0;
445 : : Assert( C(e,c,0) > -eps && C(e,c,0) < 1.0+eps,
446 : : "0 <= AEC <= 1.0 failed: C = " + std::to_string(C(e,c,0)) );
447 : : }
448 : : }
449 : :
450 : : // System limiting
451 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
452 [ + + ]: 9315492 : for (const auto& sys : m_sys) {
453 : : tk::real cs = 1.0;
454 [ + + ][ + + ]: 363840 : for (auto i : sys) if (C(e,i,0) < cs) cs = C(e,i,0);
455 [ + + ]: 363840 : for (auto i : sys) C(e,i,0) = cs;
456 : : }
457 : : }
458 : :
459 : : // save the limited antidiffusive element contributions (Lohner: AEC^c)
460 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
461 [ + - ]: 9254852 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
462 : 9254852 : inpoel[e*4+2], inpoel[e*4+3] }};
463 : :
464 : : // access pointer to solution at element nodes
465 [ + - ][ - - ]: 9254852 : std::vector< const tk::real* > a( ncomp );
466 [ + + ]: 23143444 : for (ncomp_t c=0; c<ncomp; ++c) a[c] = A.cptr( c, 0 );
467 : :
468 : : // Scatter-add limited antidiffusive element contributions to nodes. At
469 : : // nodes where Dirichlet boundary conditions are set, the AECs are set to
470 : : // zero so the limit coefficient has no effect. This yields no increment for
471 : : // those nodes. See the detailed discussion when computing the AECs.
472 [ + + ]: 46274260 : for (std::size_t j=0; j<4; ++j) {
473 : : auto b = bcdir.find( N[j] ); // Dirichlet BC
474 [ + + ]: 92573776 : for (ncomp_t c=0; c<ncomp; ++c) {
475 [ + + ][ + - ]: 55554368 : if (b != end(bcdir) && b->second[c].first) {
476 : 14546180 : A.var(a[c],N[j]) += m_aec(e*4+j,c,0);
477 : : } else {
478 : 41008188 : A.var(a[c],N[j]) += C(e,c,0) * m_aec(e*4+j,c,0);
479 : : }
480 : : }
481 : : }
482 : : }
483 : 18283 : }
484 : :
485 : : std::tuple< std::vector< std::string >,
486 : : std::vector< std::vector< tk::real > > >
487 : 2233 : FluxCorrector::fields( const std::vector< std::size_t >& /*inpoel*/ ) const
488 : : // *****************************************************************************
489 : : // Collect mesh output fields from FCT
490 : : //! \return Names and fields in mesh cells
491 : : // *****************************************************************************
492 : : {
493 : : using tuple_t = std::tuple< std::vector< std::string >,
494 : : std::vector< std::vector< tk::real > > >;
495 : 2233 : return tuple_t{};
496 : : }
|
