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 : 16320 : tk::srcInt( ncomp_t system,
24 : : ncomp_t offset,
25 : : real t,
26 : : const std::size_t ndof,
27 : : const std::size_t nelem,
28 : : const std::vector< std::size_t >& inpoel,
29 : : const UnsMesh::Coords& coord,
30 : : const Fields& geoElem,
31 : : const CompFlowSrcFn& src,
32 : : const std::vector< std::size_t >& ndofel,
33 : : Fields& R )
34 : : // *****************************************************************************
35 : : // Compute source term integrals for DG
36 : : //! \param[in] system Equation system index
37 : : //! \param[in] offset Offset this PDE system operates from
38 : : //! \param[in] t Physical time
39 : : //! \param[in] ndof Maximum number of degrees of freedom
40 : : //! \param[in] nelem Maximum number of elements
41 : : //! \param[in] inpoel Element-node connectivity
42 : : //! \param[in] coord Array of nodal coordinates
43 : : //! \param[in] geoElem Element geometry array
44 : : //! \param[in] src Source function to use
45 : : //! \param[in] ndofel Vector of local number of degrees of freedome
46 : : //! \param[in,out] R Right-hand side vector computed
47 : : // *****************************************************************************
48 : : {
49 : 16320 : const auto& cx = coord[0];
50 : 16320 : const auto& cy = coord[1];
51 : 16320 : const auto& cz = coord[2];
52 : :
53 [ + + ]: 5573700 : for (std::size_t e=0; e<nelem; ++e)
54 : : {
55 [ + - ]: 5557380 : auto ng = tk::NGvol(ndofel[e]);
56 : :
57 : : // arrays for quadrature points
58 : 11114760 : std::array< std::vector< real >, 3 > coordgp;
59 : 11114760 : std::vector< real > wgp;
60 : :
61 [ + - ]: 5557380 : coordgp[0].resize( ng );
62 [ + - ]: 5557380 : coordgp[1].resize( ng );
63 [ + - ]: 5557380 : coordgp[2].resize( ng );
64 [ + - ]: 5557380 : wgp.resize( ng );
65 : :
66 [ + - ]: 5557380 : GaussQuadratureTet( ng, coordgp, wgp );
67 : :
68 : : // Extract the element coordinates
69 : : std::array< std::array< real, 3>, 4 > coordel {{
70 : 16672140 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
71 : 16672140 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
72 : 16672140 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
73 : 50016420 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
74 : :
75 [ + + ]: 27115032 : for (std::size_t igp=0; igp<ng; ++igp)
76 : : {
77 : : // Compute the coordinates of quadrature point at physical domain
78 [ + - ]: 21557652 : auto gp = eval_gp( igp, coordel, coordgp );
79 : :
80 : : // Compute the basis function
81 : : auto B =
82 [ + - ]: 43115304 : eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
83 : :
84 : : // Compute the source term variable
85 : : std::array< real, 5 > s;
86 [ + - ]: 21557652 : src( system, gp[0], gp[1], gp[2], t, s[0], s[1], s[2], s[3], s[4] );
87 : :
88 [ + - ]: 21557652 : auto wt = wgp[igp] * geoElem(e, 0, 0);
89 : :
90 [ + - ]: 21557652 : update_rhs( offset, ndof, ndofel[e], wt, e, B, s, R );
91 : : }
92 : : }
93 : 16320 : }
94 : :
95 : : void
96 : 21557652 : tk::update_rhs( ncomp_t offset,
97 : : const std::size_t ndof,
98 : : const std::size_t ndof_el,
99 : : const tk::real wt,
100 : : const std::size_t e,
101 : : const std::vector< tk::real >& B,
102 : : const std::array< tk::real, 5 >& s,
103 : : Fields& R )
104 : : // *****************************************************************************
105 : : // Update the rhs by adding the source term integrals
106 : : //! \param[in] offset Offset this PDE system operates from
107 : : //! \param[in] ndof Maximum number of degrees of freedom
108 : : //! \param[in] ndof_el Number of degrees of freedom for local element
109 : : //! \param[in] wt Weight of gauss quadrature point
110 : : //! \param[in] e Element index
111 : : //! \param[in] B Vector of basis functions
112 : : //! \param[in] s Source terms for each of the 5 components
113 : : //! \param[in,out] R Right-hand side vector computed
114 : : // *****************************************************************************
115 : : {
116 [ - + ][ - - ]: 21557652 : Assert( B.size() == ndof_el, "Size mismatch for basis function" );
[ - - ][ - - ]
117 : :
118 [ + + ]: 129345912 : for (ncomp_t c=0; c<5; ++c)
119 : : {
120 : 107788260 : auto mark = c*ndof;
121 : 107788260 : R(e, mark, offset) += wt * s[c];
122 : :
123 [ + + ]: 107788260 : if ( ndof_el > 1 )
124 : : {
125 : 91335375 : R(e, mark+1, offset) += wt * s[c] * B[1];
126 : 91335375 : R(e, mark+2, offset) += wt * s[c] * B[2];
127 : 91335375 : R(e, mark+3, offset) += wt * s[c] * B[3];
128 : :
129 [ + + ]: 91335375 : if( ndof_el > 4 )
130 : : {
131 : 63553050 : R(e, mark+4, offset) += wt * s[c] * B[4];
132 : 63553050 : R(e, mark+5, offset) += wt * s[c] * B[5];
133 : 63553050 : R(e, mark+6, offset) += wt * s[c] * B[6];
134 : 63553050 : R(e, mark+7, offset) += wt * s[c] * B[7];
135 : 63553050 : R(e, mark+8, offset) += wt * s[c] * B[8];
136 : 63553050 : R(e, mark+9, offset) += wt * s[c] * B[9];
137 : : }
138 : : }
139 : : }
140 : 21557652 : }
|
