Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Source.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 Functions for computing integrals of an arbitrary source term of a
9 : : for single-material compressible flow, CompFlow using DG methods
10 : : \details This file contains functionality for computing integrals of an
11 : : arbitrary source term for single-material compressible flow, CompFlow with
12 : : discontinuous Galerkin methods for various orders of numerical
13 : : representation.
14 : : */
15 : : // *****************************************************************************
16 : :
17 : : #include <vector>
18 : :
19 : : #include "Source.hpp"
20 : : #include "Quadrature.hpp"
21 : :
22 : : void
23 : 20925 : tk::srcInt( const std::vector< inciter::EOS >& mat_blk,
24 : : real t,
25 : : const std::size_t ndof,
26 : : const std::size_t nelem,
27 : : const std::vector< std::size_t >& inpoel,
28 : : const UnsMesh::Coords& coord,
29 : : const Fields& geoElem,
30 : : const SrcFn& src,
31 : : const std::vector< std::size_t >& ndofel,
32 : : Fields& R,
33 : : std::size_t nmat )
34 : : // *****************************************************************************
35 : : // Compute source term integrals for DG
36 : : //! \param[in] mat_blk Material block EOS
37 : : //! \param[in] t Physical time
38 : : //! \param[in] ndof Maximum number of degrees of freedom
39 : : //! \param[in] nelem Maximum number of elements
40 : : //! \param[in] inpoel Element-node connectivity
41 : : //! \param[in] coord Array of nodal coordinates
42 : : //! \param[in] geoElem Element geometry array
43 : : //! \param[in] src Source function to use
44 : : //! \param[in] ndofel Vector of local number of degrees of freedome
45 : : //! \param[in,out] R Right-hand side vector computed
46 : : //! \param[in] nmat Number of materials. A default is set to 1, so that calling
47 : : //! code for single material systems primitive quantities does not need to
48 : : //! specify this argument.
49 : : // *****************************************************************************
50 : : {
51 : 20925 : auto ncomp = R.nprop()/ndof;
52 : 20925 : const auto& cx = coord[0];
53 : 20925 : const auto& cy = coord[1];
54 : 20925 : const auto& cz = coord[2];
55 : :
56 [ + + ]: 6848025 : for (std::size_t e=0; e<nelem; ++e)
57 : : {
58 [ + - ]: 6827100 : auto ng = tk::NGvol(ndofel[e]);
59 : :
60 : : // arrays for quadrature points
61 : 13654200 : std::array< std::vector< real >, 3 > coordgp;
62 : 13654200 : std::vector< real > wgp;
63 : :
64 [ + - ]: 6827100 : coordgp[0].resize( ng );
65 [ + - ]: 6827100 : coordgp[1].resize( ng );
66 [ + - ]: 6827100 : coordgp[2].resize( ng );
67 [ + - ]: 6827100 : wgp.resize( ng );
68 : :
69 [ + - ]: 6827100 : GaussQuadratureTet( ng, coordgp, wgp );
70 : :
71 : : // Extract the element coordinates
72 : : std::array< std::array< real, 3>, 4 > coordel {{
73 : 20481300 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
74 : 20481300 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
75 : 20481300 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
76 : 61443900 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
77 : :
78 [ + + ]: 28048464 : for (std::size_t igp=0; igp<ng; ++igp)
79 : : {
80 : : // Compute the coordinates of quadrature point at physical domain
81 [ + - ]: 21221364 : auto gp = eval_gp( igp, coordel, coordgp );
82 : :
83 : : // Compute the basis function
84 : : auto B =
85 [ + - ]: 42442728 : eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
86 : :
87 : : // Compute the source term variable
88 [ + - ]: 42442728 : std::vector< real > s(ncomp, 0.0);
89 [ + - ]: 21221364 : src( nmat, mat_blk, gp[0], gp[1], gp[2], t, s );
90 : :
91 [ + - ]: 21221364 : auto wt = wgp[igp] * geoElem(e, 0);
92 : :
93 [ + - ]: 21221364 : update_rhs( ndof, ndofel[e], wt, e, B, s, R );
94 : : }
95 : : }
96 : 20925 : }
97 : :
98 : : void
99 : 21221364 : tk::update_rhs( const std::size_t ndof,
100 : : const std::size_t ndof_el,
101 : : const tk::real wt,
102 : : const std::size_t e,
103 : : const std::vector< tk::real >& B,
104 : : const std::vector< tk::real >& s,
105 : : Fields& R )
106 : : // *****************************************************************************
107 : : // Update the rhs by adding the source term integrals
108 : : //! \param[in] ndof Maximum number of degrees of freedom
109 : : //! \param[in] ndof_el Number of degrees of freedom for local element
110 : : //! \param[in] wt Weight of gauss quadrature point
111 : : //! \param[in] e Element index
112 : : //! \param[in] B Vector of basis functions
113 : : //! \param[in] s Source term vector
114 : : //! \param[in,out] R Right-hand side vector computed
115 : : // *****************************************************************************
116 : : {
117 [ - + ][ - - ]: 21221364 : Assert( B.size() == ndof_el, "Size mismatch for basis function" );
[ - - ][ - - ]
118 : :
119 [ + + ]: 147258084 : for (ncomp_t c=0; c<s.size(); ++c)
120 : : {
121 : 126036720 : auto mark = c*ndof;
122 : 126036720 : R(e, mark) += wt * s[c];
123 : :
124 [ + + ]: 126036720 : if ( ndof_el > 1 )
125 : : {
126 : 94040505 : R(e, mark+1) += wt * s[c] * B[1];
127 : 94040505 : R(e, mark+2) += wt * s[c] * B[2];
128 : 94040505 : R(e, mark+3) += wt * s[c] * B[3];
129 : :
130 [ + + ]: 94040505 : if( ndof_el > 4 )
131 : : {
132 : 43481130 : R(e, mark+4) += wt * s[c] * B[4];
133 : 43481130 : R(e, mark+5) += wt * s[c] * B[5];
134 : 43481130 : R(e, mark+6) += wt * s[c] * B[6];
135 : 43481130 : R(e, mark+7) += wt * s[c] * B[7];
136 : 43481130 : R(e, mark+8) += wt * s[c] * B[8];
137 : 43481130 : R(e, mark+9) += wt * s[c] * B[9];
138 : : }
139 : : }
140 : : }
141 : 21221364 : }
142 : :
143 : : void
144 : 1668 : tk::srcIntFV( const std::vector< inciter::EOS >& mat_blk,
145 : : real t,
146 : : const std::size_t nelem,
147 : : const Fields& geoElem,
148 : : const SrcFn& src,
149 : : Fields& R,
150 : : std::size_t nmat )
151 : : // *****************************************************************************
152 : : // Compute source term integrals for DG
153 : : //! \param[in] mat_blk Material block EOS
154 : : //! \param[in] t Physical time
155 : : //! \param[in] nelem Maximum number of elements
156 : : //! \param[in] geoElem Element geometry array
157 : : //! \param[in] src Source function to use
158 : : //! \param[in,out] R Right-hand side vector computed
159 : : //! \param[in] nmat Number of materials. A default is set to 1, so that calling
160 : : //! code for single material systems primitive quantities does not need to
161 : : //! specify this argument.
162 : : // *****************************************************************************
163 : : {
164 : 1668 : auto ncomp = R.nprop();
165 : :
166 [ + + ]: 298228 : for (std::size_t e=0; e<nelem; ++e)
167 : : {
168 : : // Compute the source term variable
169 [ + - ]: 593120 : std::vector< real > s(ncomp, 0.0);
170 [ + - ][ + - ]: 296560 : src( nmat, mat_blk, geoElem(e,1), geoElem(e,2), geoElem(e,3), t, s );
[ + - ][ + - ]
171 : :
172 : : // Add the source term to the rhs
173 [ + + ]: 3123280 : for (ncomp_t c=0; c<ncomp; ++c)
174 : : {
175 [ + - ][ + - ]: 2826720 : R(e, c) += geoElem(e,0) * s[c];
176 : : }
177 : : }
178 : 1668 : }
|
