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 : 34335 : tk::volInt( ncomp_t system,
23 : : std::size_t nmat,
24 : : ncomp_t offset,
25 : : real t,
26 : : const std::size_t ndof,
27 : : const std::size_t rdof,
28 : : const std::size_t nelem,
29 : : const std::vector< std::size_t >& inpoel,
30 : : const UnsMesh::Coords& coord,
31 : : const Fields& geoElem,
32 : : const FluxFn& flux,
33 : : const VelFn& vel,
34 : : const Fields& U,
35 : : const Fields& P,
36 : : const std::vector< std::size_t >& ndofel,
37 : : Fields& R,
38 : : int intsharp )
39 : : // *****************************************************************************
40 : : // Compute volume integrals for DG
41 : : //! \param[in] system Equation system index
42 : : //! \param[in] nmat Number of materials in this PDE system
43 : : //! \param[in] offset Offset this PDE system operates from
44 : : //! \param[in] t Physical time
45 : : //! \param[in] ndof Maximum number of degrees of freedom
46 : : //! \param[in] rdof Total number of degrees of freedom included reconstructed ones
47 : : //! \param[in] nelem Maximum number of elements
48 : : //! \param[in] inpoel Element-node connectivity
49 : : //! \param[in] coord Array of nodal coordinates
50 : : //! \param[in] geoElem Element geometry array
51 : : //! \param[in] flux Flux function to use
52 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
53 : : //! \param[in] U Solution vector at recent time step
54 : : //! \param[in] P Vector of primitives at recent time step
55 : : //! \param[in] ndofel Vector of local number of degrees of freedom
56 : : //! \param[in,out] R Right-hand side vector added to
57 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
58 : : //! default 0, so that it is unused for single-material and transport.
59 : : // *****************************************************************************
60 : : {
61 : 34335 : const auto& cx = coord[0];
62 : 34335 : const auto& cy = coord[1];
63 : 34335 : const auto& cz = coord[2];
64 : :
65 : 34335 : auto ncomp = U.nprop()/rdof;
66 : 34335 : auto nprim = P.nprop()/rdof;
67 : :
68 : : // compute volume integrals
69 [ + + ]: 6953415 : for (std::size_t e=0; e<nelem; ++e)
70 : : {
71 [ + + ]: 6919080 : if(ndofel[e] > 1)
72 : : {
73 [ + - ]: 6000303 : auto ng = tk::NGvol(ndofel[e]);
74 : :
75 : : // arrays for quadrature points
76 : 12000606 : std::array< std::vector< real >, 3 > coordgp;
77 : 12000606 : std::vector< real > wgp;
78 : :
79 [ + - ]: 6000303 : coordgp[0].resize( ng );
80 [ + - ]: 6000303 : coordgp[1].resize( ng );
81 [ + - ]: 6000303 : coordgp[2].resize( ng );
82 [ + - ]: 6000303 : wgp.resize( ng );
83 : :
84 [ + - ]: 6000303 : GaussQuadratureTet( ng, coordgp, wgp );
85 : :
86 : : // Extract the element coordinates
87 : : std::array< std::array< real, 3>, 4 > coordel {{
88 : 18000909 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
89 : 18000909 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
90 : 18000909 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
91 : 18000909 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
92 : 66003333 : }};
93 : :
94 : : auto jacInv =
95 : 6000303 : inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
96 : :
97 : : // Compute the derivatives of basis function for DG(P1)
98 [ + - ]: 12000606 : auto dBdx = eval_dBdx_p1( ndofel[e], jacInv );
99 : :
100 : : // Gaussian quadrature
101 [ + + ]: 46213578 : for (std::size_t igp=0; igp<ng; ++igp)
102 : : {
103 [ + + ]: 40213275 : if (ndofel[e] > 4)
104 [ + - ]: 18721560 : eval_dBdx_p2( igp, coordgp, jacInv, dBdx );
105 : :
106 : : // Compute the coordinates of quadrature point at physical domain
107 [ + - ]: 40213275 : auto gp = eval_gp( igp, coordel, coordgp );
108 : :
109 : : // Compute the basis function
110 : 80426550 : auto B = eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp],
111 [ + - ]: 160853100 : coordgp[2][igp] );
112 : :
113 [ + - ]: 40213275 : auto wt = wgp[igp] * geoElem(e, 0, 0);
114 : :
115 : : auto state = evalPolynomialSol(system, offset, intsharp, ncomp, nprim,
116 : 40213275 : rdof, nmat, e, ndofel[e], inpoel, coord, geoElem,
117 [ + - ]: 80426550 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
118 : :
119 : : // evaluate prescribed velocity (if any)
120 [ + - ]: 80426550 : auto v = vel( system, ncomp, gp[0], gp[1], gp[2], t );
121 : :
122 : : // comput flux
123 [ + - ]: 80426550 : auto fl = flux( system, ncomp, state, v );
124 : :
125 [ + - ]: 40213275 : update_rhs( ncomp, offset, ndof, ndofel[e], wt, e, dBdx, fl, R );
126 : : }
127 : : }
128 : : }
129 : 34335 : }
130 : :
131 : : void
132 : 40213275 : tk::update_rhs( ncomp_t ncomp,
133 : : ncomp_t offset,
134 : : const std::size_t ndof,
135 : : const std::size_t ndof_el,
136 : : const tk::real wt,
137 : : const std::size_t e,
138 : : const std::array< std::vector<tk::real>, 3 >& dBdx,
139 : : const std::vector< std::array< tk::real, 3 > >& fl,
140 : : Fields& R )
141 : : // *****************************************************************************
142 : : // Update the rhs by adding the source term integrals
143 : : //! \param[in] ncomp Number of scalar components in this PDE system
144 : : //! \param[in] offset Offset this PDE system operates from
145 : : //! \param[in] ndof Maximum number of degrees of freedom
146 : : //! \param[in] ndof_el Number of degrees of freedom for local element
147 : : //! \param[in] wt Weight of gauss quadrature point
148 : : //! \param[in] e Element index
149 : : //! \param[in] dBdx Vector of basis function derivatives
150 : : //! \param[in] fl Vector of numerical flux
151 : : //! \param[in,out] R Right-hand side vector computed
152 : : // *****************************************************************************
153 : : {
154 [ - + ][ - - ]: 40213275 : Assert( dBdx[0].size() == ndof_el,
[ - - ][ - - ]
155 : : "Size mismatch for basis function derivatives" );
156 [ - + ][ - - ]: 40213275 : Assert( dBdx[1].size() == ndof_el,
[ - - ][ - - ]
157 : : "Size mismatch for basis function derivatives" );
158 [ - + ][ - - ]: 40213275 : Assert( dBdx[2].size() == ndof_el,
[ - - ][ - - ]
159 : : "Size mismatch for basis function derivatives" );
160 [ - + ][ - - ]: 40213275 : Assert( fl.size() == ncomp, "Size mismatch for flux term" );
[ - - ][ - - ]
161 : :
162 [ + + ]: 171686850 : for (ncomp_t c=0; c<ncomp; ++c)
163 : : {
164 : 131473575 : auto mark = c*ndof;
165 : 262947150 : R(e, mark+1, offset) +=
166 : 131473575 : wt * (fl[c][0]*dBdx[0][1] + fl[c][1]*dBdx[1][1] + fl[c][2]*dBdx[2][1]);
167 : 262947150 : R(e, mark+2, offset) +=
168 : 131473575 : wt * (fl[c][0]*dBdx[0][2] + fl[c][1]*dBdx[1][2] + fl[c][2]*dBdx[2][2]);
169 : 262947150 : R(e, mark+3, offset) +=
170 : 131473575 : wt * (fl[c][0]*dBdx[0][3] + fl[c][1]*dBdx[1][3] + fl[c][2]*dBdx[2][3]);
171 : :
172 [ + + ]: 131473575 : if( ndof_el > 4 )
173 : : {
174 : 139128000 : R(e, mark+4, offset) +=
175 : 69564000 : wt * (fl[c][0]*dBdx[0][4] + fl[c][1]*dBdx[1][4] + fl[c][2]*dBdx[2][4]);
176 : 139128000 : R(e, mark+5, offset) +=
177 : 69564000 : wt * (fl[c][0]*dBdx[0][5] + fl[c][1]*dBdx[1][5] + fl[c][2]*dBdx[2][5]);
178 : 139128000 : R(e, mark+6, offset) +=
179 : 69564000 : wt * (fl[c][0]*dBdx[0][6] + fl[c][1]*dBdx[1][6] + fl[c][2]*dBdx[2][6]);
180 : 139128000 : R(e, mark+7, offset) +=
181 : 69564000 : wt * (fl[c][0]*dBdx[0][7] + fl[c][1]*dBdx[1][7] + fl[c][2]*dBdx[2][7]);
182 : 139128000 : R(e, mark+8, offset) +=
183 : 69564000 : wt * (fl[c][0]*dBdx[0][8] + fl[c][1]*dBdx[1][8] + fl[c][2]*dBdx[2][8]);
184 : 69564000 : R(e, mark+9, offset) +=
185 : 69564000 : wt * (fl[c][0]*dBdx[0][9] + fl[c][1]*dBdx[1][9] + fl[c][2]*dBdx[2][9]);
186 : : }
187 : : }
188 : 40213275 : }
|
