Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Initialize.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 initialization of system of PDEs in DG methods
9 : : \details This file contains functionality for setting initial conditions
10 : : and evaluating known solutions used in discontinuous Galerkin methods for
11 : : various orders of numerical representation.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <array>
16 : : #include <vector>
17 : : #include <iostream>
18 : : #include "Data.hpp"
19 : : #include "Initialize.hpp"
20 : : #include "Quadrature.hpp"
21 : : #include "Inciter/InputDeck/InputDeck.hpp"
22 : :
23 : : namespace inciter {
24 : :
25 : : extern ctr::InputDeck g_inputdeck;
26 : :
27 : : } // inciter::
28 : :
29 : : void
30 : 1231 : tk::initialize( ncomp_t ncomp,
31 : : const std::vector< inciter::EOS >& mat_blk,
32 : : const Fields& L,
33 : : const std::vector< std::size_t >& inpoel,
34 : : const UnsMesh::Coords& coord,
35 : : const InitializeFn& solution,
36 : : Fields& unk,
37 : : real t,
38 : : const std::size_t nielem )
39 : : // *****************************************************************************
40 : : //! Initalize a system of DGPDEs by projecting the exact solution in the DG
41 : : //! solution space
42 : : //! \details This is the public interface exposed to client code.
43 : : //! \param[in] ncomp Number of scalar components in this PDE system
44 : : //! \param[in] L Block diagonal mass matrix
45 : : //! \param[in] inpoel Element-node connectivity
46 : : //! \param[in] coord Array of node coordinates
47 : : //! \param[in] solution Function to call to evaluate known solution or initial
48 : : //! conditions at x,y,z,t
49 : : //! \param[in,out] unk Array of unknowns
50 : : //! \param[in] t Physical time
51 : : //! \param[in] nielem Number of internal elements
52 : : // *****************************************************************************
53 : : {
54 : 1231 : const auto ndof = inciter::g_inputdeck.get< tag::ndof >();
55 : 1231 : const auto rdof = inciter::g_inputdeck.get< tag::rdof >();
56 : :
57 : : // Number of quadrature points for volume integration
58 : 1231 : auto ng = tk::NGinit(ndof);
59 : :
60 : : // arrays for quadrature points
61 : : std::array< std::vector< real >, 3 > coordgp;
62 : : std::vector< real > wgp;
63 : :
64 [ + - ]: 1231 : coordgp[0].resize( ng );
65 [ + - ]: 1231 : coordgp[1].resize( ng );
66 [ + - ]: 1231 : coordgp[2].resize( ng );
67 [ + - ]: 1231 : wgp.resize( ng );
68 : :
69 : : // get quadrature point weights and coordinates for triangle
70 [ + - ]: 1231 : GaussQuadratureTet( ng, coordgp, wgp );
71 : :
72 : : const auto& cx = coord[0];
73 : : const auto& cy = coord[1];
74 : : const auto& cz = coord[2];
75 : :
76 [ + + ]: 524337 : for (std::size_t e=0; e<nielem; ++e) { // for all tets
77 : : // The volume of tetrahedron
78 : 523106 : auto vole = L(e, 0);
79 : :
80 : : // Extract the element coordinates
81 : : std::array< std::array< real, 3>, 4 > coordel {{
82 : 523106 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
83 : 523106 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
84 : 523106 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
85 : 523106 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
86 : :
87 : : // right hand side vector
88 [ + - ][ - - ]: 523106 : std::vector< real > R( ncomp*ndof, 0.0 );
89 : :
90 : : // Gaussian quadrature
91 [ + + ]: 2133077 : for (std::size_t igp=0; igp<ng; ++igp)
92 : : {
93 : : // Compute the coordinates of quadrature point at physical domain
94 [ + - ]: 1609971 : auto gp = eval_gp( igp, coordel, coordgp );
95 : :
96 : : // Compute the basis function
97 : : auto B =
98 [ + - ]: 1609971 : eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
99 : :
100 [ - + ]: 1609971 : const auto s = solution( ncomp, mat_blk, gp[0], gp[1], gp[2], t );
101 : :
102 [ + - ]: 1609971 : auto wt = wgp[igp] * vole;
103 : :
104 [ + - ]: 1609971 : update_rhs( ncomp, ndof, wt, B, s, R );
105 : : }
106 : :
107 : : // Compute the initial conditions
108 [ + - ]: 523106 : eval_init(ncomp, ndof, rdof, e, R, L, unk);
109 : : }
110 : 1231 : }
111 : :
112 : : void
113 : 1609971 : tk::update_rhs( ncomp_t ncomp,
114 : : const std::size_t ndof,
115 : : const tk::real wt,
116 : : const std::vector< tk::real >& B,
117 : : const std::vector< tk::real >& s,
118 : : std::vector< tk::real >& R )
119 : : // *****************************************************************************
120 : : // Update the rhs by adding the initial analytical solution term
121 : : //! \param[in] ncomp Number of scalar components in this PDE system
122 : : //! \param[in] ndof Number of degrees of freedom
123 : : //! \param[in] wt Weight of gauss quadrature point
124 : : //! \param[in] B Vector of basis functions
125 : : //! \param[in] s Vector of analytical solution at quadrature point
126 : : //! \param[in,out] R Right-hand side vector
127 : : // *****************************************************************************
128 : : {
129 : : Assert( B.size() == ndof, "Size mismatch for basis function" );
130 : : Assert( s.size() >= ncomp, "Size mismatch for source term" );
131 : :
132 [ + + ]: 7387948 : for (ncomp_t c=0; c<ncomp; ++c)
133 : : {
134 : : // DG(P0)
135 : 5777977 : auto mark = c*ndof;
136 [ + + ]: 5777977 : R[mark] += wt * s[c];
137 : :
138 [ + + ]: 5777977 : if(ndof > 1) //DG(P1)
139 : : {
140 [ + + ]: 4848774 : R[mark+1] += wt * s[c] * B[1];
141 : 4848774 : R[mark+2] += wt * s[c] * B[2];
142 : 4848774 : R[mark+3] += wt * s[c] * B[3];
143 : :
144 [ + + ]: 4848774 : if(ndof > 4) //DG(P2)
145 : : {
146 : 2873164 : R[mark+4] += wt * s[c] * B[4];
147 : 2873164 : R[mark+5] += wt * s[c] * B[5];
148 : 2873164 : R[mark+6] += wt * s[c] * B[6];
149 : 2873164 : R[mark+7] += wt * s[c] * B[7];
150 : 2873164 : R[mark+8] += wt * s[c] * B[8];
151 : 2873164 : R[mark+9] += wt * s[c] * B[9];
152 : : }
153 : : }
154 : : }
155 : 1609971 : }
156 : :
157 : : void
158 : 523106 : tk::eval_init( ncomp_t ncomp,
159 : : const std::size_t ndof,
160 : : const std::size_t rdof,
161 : : const std::size_t e,
162 : : const std::vector< tk::real >& R,
163 : : const Fields& L,
164 : : Fields& unk )
165 : : // *****************************************************************************
166 : : // Compute the initial conditions
167 : : //! \param[in] ncomp Number of scalar components in this PDE system
168 : : //! \param[in] ndof Number of degrees of freedom
169 : : //! \param[in] rdof Total number of reconstructed degrees of freedom
170 : : //! \param[in] e Element index
171 : : //! \param[in] R Right-hand side vector
172 : : //! \param[in] L Block diagonal mass matrix
173 : : //! \param[in,out] unk Array of unknowns
174 : : // *****************************************************************************
175 : : {
176 [ + + ]: 1798650 : for (ncomp_t c=0; c<ncomp; ++c)
177 : : {
178 : : // DG(P0)
179 : 1275544 : auto mark = c*ndof;
180 : 1275544 : auto rmark = c*rdof;
181 [ + + ]: 1275544 : unk(e, rmark) = R[mark] / L(e, mark);
182 : :
183 : : // if P0P1, initialize higher dofs to 0
184 [ + + ]: 1275544 : if (rdof > ndof)
185 : : {
186 : 234216 : unk(e, rmark+1) = 0.0;
187 : 234216 : unk(e, rmark+2) = 0.0;
188 : 234216 : unk(e, rmark+3) = 0.0;
189 : : }
190 : :
191 [ + + ]: 1275544 : if(ndof > 1) // DG(P1)
192 : : {
193 [ + + ]: 346341 : unk(e, rmark+1) = R[mark+1] / L(e, mark+1);
194 : 346341 : unk(e, rmark+2) = R[mark+2] / L(e, mark+2);
195 : 346341 : unk(e, rmark+3) = R[mark+3] / L(e, mark+3);
196 : :
197 [ + + ]: 346341 : if(ndof > 4) // DG(P2)
198 : : {
199 : 205226 : unk(e, rmark+4) = R[mark+4] / L(e, mark+4);
200 : 205226 : unk(e, rmark+5) = R[mark+5] / L(e, mark+5);
201 : 205226 : unk(e, rmark+6) = R[mark+6] / L(e, mark+6);
202 : 205226 : unk(e, rmark+7) = R[mark+7] / L(e, mark+7);
203 : 205226 : unk(e, rmark+8) = R[mark+8] / L(e, mark+8);
204 : 205226 : unk(e, rmark+9) = R[mark+9] / L(e, mark+9);
205 : : }
206 : : }
207 : : }
208 : 523106 : }
|
