Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Volume.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 volume integrals for a system of PDEs in DG
9 : : methods
10 : : \details This file contains functionality for computing volume integrals for
11 : : a system of PDEs used in discontinuous Galerkin methods for various orders
12 : : of numerical representation.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include "Volume.hpp"
17 : : #include "Vector.hpp"
18 : : #include "Quadrature.hpp"
19 : : #include "Reconstruction.hpp"
20 : :
21 : : void
22 : 34065 : tk::volInt( std::size_t nmat,
23 : : real t,
24 : : const std::vector< inciter::EOS >& mat_blk,
25 : : const std::size_t ndof,
26 : : const std::size_t rdof,
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 FluxFn& flux,
32 : : const VelFn& vel,
33 : : const Fields& U,
34 : : const Fields& P,
35 : : const std::vector< std::size_t >& ndofel,
36 : : Fields& R,
37 : : int intsharp )
38 : : // *****************************************************************************
39 : : // Compute volume integrals for DG
40 : : //! \param[in] nmat Number of materials in this PDE system
41 : : //! \param[in] t Physical time
42 : : //! \param[in] mat_blk EOS material block
43 : : //! \param[in] ndof Maximum number of degrees of freedom
44 : : //! \param[in] rdof Total number of degrees of freedom included reconstructed ones
45 : : //! \param[in] nelem Maximum number of elements
46 : : //! \param[in] inpoel Element-node connectivity
47 : : //! \param[in] coord Array of nodal coordinates
48 : : //! \param[in] geoElem Element geometry array
49 : : //! \param[in] flux Flux function to use
50 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
51 : : //! \param[in] U Solution vector at recent time step
52 : : //! \param[in] P Vector of primitives at recent time step
53 : : //! \param[in] ndofel Vector of local number of degrees of freedom
54 : : //! \param[in,out] R Right-hand side vector added to
55 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
56 : : //! default 0, so that it is unused for single-material and transport.
57 : : // *****************************************************************************
58 : : {
59 : 34065 : const auto& cx = coord[0];
60 : 34065 : const auto& cy = coord[1];
61 : 34065 : const auto& cz = coord[2];
62 : :
63 : 34065 : auto ncomp = U.nprop()/rdof;
64 : 34065 : auto nprim = P.nprop()/rdof;
65 : :
66 : : // compute volume integrals
67 [ + + ]: 6561465 : for (std::size_t e=0; e<nelem; ++e)
68 : : {
69 [ + + ]: 6527400 : if(ndofel[e] > 1)
70 : : {
71 [ + - ]: 5691417 : auto ng = tk::NGvol(ndofel[e]);
72 : :
73 : : // arrays for quadrature points
74 : 11382834 : std::array< std::vector< real >, 3 > coordgp;
75 : 11382834 : std::vector< real > wgp;
76 : :
77 [ + - ]: 5691417 : coordgp[0].resize( ng );
78 [ + - ]: 5691417 : coordgp[1].resize( ng );
79 [ + - ]: 5691417 : coordgp[2].resize( ng );
80 [ + - ]: 5691417 : wgp.resize( ng );
81 : :
82 [ + - ]: 5691417 : GaussQuadratureTet( ng, coordgp, wgp );
83 : :
84 : : // Extract the element coordinates
85 : : std::array< std::array< real, 3>, 4 > coordel {{
86 : 17074251 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
87 : 17074251 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
88 : 17074251 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
89 : 17074251 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
90 : 62605587 : }};
91 : :
92 : : auto jacInv =
93 : 5691417 : inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
94 : :
95 : : // Compute the derivatives of basis function for DG(P1)
96 [ + - ]: 11382834 : auto dBdx = eval_dBdx_p1( ndofel[e], jacInv );
97 : :
98 : : // Gaussian quadrature
99 [ + + ]: 42170598 : for (std::size_t igp=0; igp<ng; ++igp)
100 : : {
101 [ + + ]: 36479181 : if (ndofel[e] > 4)
102 [ + - ]: 14707176 : eval_dBdx_p2( igp, coordgp, jacInv, dBdx );
103 : :
104 : : // Compute the coordinates of quadrature point at physical domain
105 [ + - ]: 36479181 : auto gp = eval_gp( igp, coordel, coordgp );
106 : :
107 : : // Compute the basis function
108 : 72958362 : auto B = eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp],
109 [ + - ]: 145916724 : coordgp[2][igp] );
110 : :
111 [ + - ]: 36479181 : auto wt = wgp[igp] * geoElem(e, 0);
112 : :
113 : : auto state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
114 : 36479181 : rdof, nmat, e, ndofel[e], inpoel, coord, geoElem,
115 [ + - ]: 72958362 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
116 : :
117 : : // evaluate prescribed velocity (if any)
118 [ + - ]: 72958362 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
119 : :
120 : : // comput flux
121 [ + - ]: 72958362 : auto fl = flux( ncomp, mat_blk, state, v );
122 : :
123 [ + - ]: 36479181 : update_rhs( ncomp, ndof, ndofel[e], wt, e, dBdx, fl, R );
124 : : }
125 : : }
126 : : }
127 : 34065 : }
128 : :
129 : : void
130 : 36479181 : tk::update_rhs( ncomp_t ncomp,
131 : : const std::size_t ndof,
132 : : const std::size_t ndof_el,
133 : : const tk::real wt,
134 : : const std::size_t e,
135 : : const std::array< std::vector<tk::real>, 3 >& dBdx,
136 : : const std::vector< std::array< tk::real, 3 > >& fl,
137 : : Fields& R )
138 : : // *****************************************************************************
139 : : // Update the rhs by adding the source term integrals
140 : : //! \param[in] ncomp Number of scalar components in this PDE system
141 : : //! \param[in] ndof Maximum number of degrees of freedom
142 : : //! \param[in] ndof_el Number of degrees of freedom for local element
143 : : //! \param[in] wt Weight of gauss quadrature point
144 : : //! \param[in] e Element index
145 : : //! \param[in] dBdx Vector of basis function derivatives
146 : : //! \param[in] fl Vector of numerical flux
147 : : //! \param[in,out] R Right-hand side vector computed
148 : : // *****************************************************************************
149 : : {
150 [ - + ][ - - ]: 36479181 : Assert( dBdx[0].size() == ndof_el,
[ - - ][ - - ]
151 : : "Size mismatch for basis function derivatives" );
152 [ - + ][ - - ]: 36479181 : Assert( dBdx[1].size() == ndof_el,
[ - - ][ - - ]
153 : : "Size mismatch for basis function derivatives" );
154 [ - + ][ - - ]: 36479181 : Assert( dBdx[2].size() == ndof_el,
[ - - ][ - - ]
155 : : "Size mismatch for basis function derivatives" );
156 [ - + ][ - - ]: 36479181 : Assert( fl.size() == ncomp, "Size mismatch for flux term" );
[ - - ][ - - ]
157 : :
158 [ + + ]: 150191886 : for (ncomp_t c=0; c<ncomp; ++c)
159 : : {
160 : 113712705 : auto mark = c*ndof;
161 : 227425410 : R(e, mark+1) +=
162 : 113712705 : wt * (fl[c][0]*dBdx[0][1] + fl[c][1]*dBdx[1][1] + fl[c][2]*dBdx[2][1]);
163 : 227425410 : R(e, mark+2) +=
164 : 113712705 : wt * (fl[c][0]*dBdx[0][2] + fl[c][1]*dBdx[1][2] + fl[c][2]*dBdx[2][2]);
165 : 227425410 : R(e, mark+3) +=
166 : 113712705 : wt * (fl[c][0]*dBdx[0][3] + fl[c][1]*dBdx[1][3] + fl[c][2]*dBdx[2][3]);
167 : :
168 [ + + ]: 113712705 : if( ndof_el > 4 )
169 : : {
170 : 98984160 : R(e, mark+4) +=
171 : 49492080 : wt * (fl[c][0]*dBdx[0][4] + fl[c][1]*dBdx[1][4] + fl[c][2]*dBdx[2][4]);
172 : 98984160 : R(e, mark+5) +=
173 : 49492080 : wt * (fl[c][0]*dBdx[0][5] + fl[c][1]*dBdx[1][5] + fl[c][2]*dBdx[2][5]);
174 : 98984160 : R(e, mark+6) +=
175 : 49492080 : wt * (fl[c][0]*dBdx[0][6] + fl[c][1]*dBdx[1][6] + fl[c][2]*dBdx[2][6]);
176 : 98984160 : R(e, mark+7) +=
177 : 49492080 : wt * (fl[c][0]*dBdx[0][7] + fl[c][1]*dBdx[1][7] + fl[c][2]*dBdx[2][7]);
178 : 98984160 : R(e, mark+8) +=
179 : 49492080 : wt * (fl[c][0]*dBdx[0][8] + fl[c][1]*dBdx[1][8] + fl[c][2]*dBdx[2][8]);
180 : 49492080 : R(e, mark+9) +=
181 : 49492080 : wt * (fl[c][0]*dBdx[0][9] + fl[c][1]*dBdx[1][9] + fl[c][2]*dBdx[2][9]);
182 : : }
183 : : }
184 : 36479181 : }
|
