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 : 1050 : tk::initialize( ncomp_t system,
31 : : ncomp_t ncomp,
32 : : ncomp_t offset,
33 : : const Fields& L,
34 : : const std::vector< std::size_t >& inpoel,
35 : : const UnsMesh::Coords& coord,
36 : : const InitializeFn& solution,
37 : : Fields& unk,
38 : : real t,
39 : : const std::size_t nielem )
40 : : // *****************************************************************************
41 : : //! Initalize a system of DGPDEs by projecting the exact solution in the DG
42 : : //! solution space
43 : : //! \details This is the public interface exposed to client code.
44 : : //! \param[in] system Equation system index
45 : : //! \param[in] ncomp Number of scalar components in this PDE system
46 : : //! \param[in] offset Offset this PDE system operates from
47 : : //! \param[in] L Block diagonal mass matrix
48 : : //! \param[in] inpoel Element-node connectivity
49 : : //! \param[in] coord Array of node coordinates
50 : : //! \param[in] solution Function to call to evaluate known solution or initial
51 : : //! conditions at x,y,z,t
52 : : //! \param[in,out] unk Array of unknowns
53 : : //! \param[in] t Physical time
54 : : //! \param[in] nielem Number of internal elements
55 : : // *****************************************************************************
56 : : {
57 : 1050 : const auto ndof = inciter::g_inputdeck.get< tag::discr, tag::ndof >();
58 : 1050 : const auto rdof = inciter::g_inputdeck.get< tag::discr, tag::rdof >();
59 : :
60 : : // Number of quadrature points for volume integration
61 [ + - ]: 1050 : auto ng = tk::NGinit(ndof);
62 : :
63 : : // arrays for quadrature points
64 : 2100 : std::array< std::vector< real >, 3 > coordgp;
65 : 2100 : std::vector< real > wgp;
66 : :
67 [ + - ]: 1050 : coordgp[0].resize( ng );
68 [ + - ]: 1050 : coordgp[1].resize( ng );
69 [ + - ]: 1050 : coordgp[2].resize( ng );
70 [ + - ]: 1050 : wgp.resize( ng );
71 : :
72 : : // get quadrature point weights and coordinates for triangle
73 [ + - ]: 1050 : GaussQuadratureTet( ng, coordgp, wgp );
74 : :
75 : 1050 : const auto& cx = coord[0];
76 : 1050 : const auto& cy = coord[1];
77 : 1050 : const auto& cz = coord[2];
78 : :
79 [ + + ]: 1078019 : for (std::size_t e=0; e<nielem; ++e) { // for all tets
80 : : // The volume of tetrahedron
81 [ + - ]: 1076969 : auto vole = L(e, 0, offset);
82 : :
83 : : // Extract the element coordinates
84 : : std::array< std::array< real, 3>, 4 > coordel {{
85 : 3230907 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
86 : 3230907 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
87 : 3230907 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
88 : 9692721 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }} }};
89 : :
90 : : // right hand side vector
91 [ + - ]: 2153938 : std::vector< real > R( ncomp*ndof, 0.0 );
92 : :
93 : : // Gaussian quadrature
94 [ + + ]: 3315813 : for (std::size_t igp=0; igp<ng; ++igp)
95 : : {
96 : : // Compute the coordinates of quadrature point at physical domain
97 [ + - ]: 2238844 : auto gp = eval_gp( igp, coordel, coordgp );
98 : :
99 : : // Compute the basis function
100 : : auto B =
101 [ + - ]: 4477688 : eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
102 : :
103 [ + - ]: 4477688 : const auto s = solution( system, ncomp, gp[0], gp[1], gp[2], t );
104 : :
105 : 2238844 : auto wt = wgp[igp] * vole;
106 : :
107 [ + - ]: 2238844 : update_rhs( ncomp, ndof, wt, B, s, R );
108 : : }
109 : :
110 : : // Compute the initial conditions
111 [ + - ]: 1076969 : eval_init(ncomp, offset, ndof, rdof, e, R, L, unk);
112 : : }
113 : 1050 : }
114 : :
115 : : void
116 : 2238844 : tk::update_rhs( ncomp_t ncomp,
117 : : const std::size_t ndof,
118 : : const tk::real wt,
119 : : const std::vector< tk::real >& B,
120 : : const std::vector< tk::real >& s,
121 : : std::vector< tk::real >& R )
122 : : // *****************************************************************************
123 : : // Update the rhs by adding the initial analytical solution term
124 : : //! \param[in] ncomp Number of scalar components in this PDE system
125 : : //! \param[in] ndof Number of degrees of freedom
126 : : //! \param[in] wt Weight of gauss quadrature point
127 : : //! \param[in] B Vector of basis functions
128 : : //! \param[in] s Vector of analytical solution at quadrature point
129 : : //! \param[in,out] R Right-hand side vector
130 : : // *****************************************************************************
131 : : {
132 [ - + ][ - - ]: 2238844 : Assert( B.size() == ndof, "Size mismatch for basis function" );
[ - - ][ - - ]
133 [ - + ][ - - ]: 2238844 : Assert( s.size() >= ncomp, "Size mismatch for source term" );
[ - - ][ - - ]
134 : :
135 [ + + ]: 8725922 : for (ncomp_t c=0; c<ncomp; ++c)
136 : : {
137 : : // DG(P0)
138 : 6487078 : auto mark = c*ndof;
139 : 6487078 : R[mark] += wt * s[c];
140 : :
141 [ + + ]: 6487078 : if(ndof > 1) //DG(P1)
142 : : {
143 : 5167778 : R[mark+1] += wt * s[c] * B[1];
144 : 5167778 : R[mark+2] += wt * s[c] * B[2];
145 : 5167778 : R[mark+3] += wt * s[c] * B[3];
146 : :
147 [ + + ]: 5167778 : if(ndof > 4) //DG(P2)
148 : : {
149 : 3383184 : R[mark+4] += wt * s[c] * B[4];
150 : 3383184 : R[mark+5] += wt * s[c] * B[5];
151 : 3383184 : R[mark+6] += wt * s[c] * B[6];
152 : 3383184 : R[mark+7] += wt * s[c] * B[7];
153 : 3383184 : R[mark+8] += wt * s[c] * B[8];
154 : 3383184 : R[mark+9] += wt * s[c] * B[9];
155 : : }
156 : : }
157 : : }
158 : 2238844 : }
159 : :
160 : : void
161 : 1076969 : tk::eval_init( ncomp_t ncomp,
162 : : ncomp_t offset,
163 : : const std::size_t ndof,
164 : : const std::size_t rdof,
165 : : const std::size_t e,
166 : : const std::vector< tk::real >& R,
167 : : const Fields& L,
168 : : Fields& unk )
169 : : // *****************************************************************************
170 : : // Compute the initial conditions
171 : : //! \param[in] ncomp Number of scalar components in this PDE system
172 : : //! \param[in] offset Offset this PDE system operates from
173 : : //! \param[in] ndof Number of degrees of freedom
174 : : //! \param[in] rdof Total number of reconstructed degrees of freedom
175 : : //! \param[in] e Element index
176 : : //! \param[in] R Right-hand side vector
177 : : //! \param[in] L Block diagonal mass matrix
178 : : //! \param[in,out] unk Array of unknowns
179 : : // *****************************************************************************
180 : : {
181 [ + + ]: 2765396 : for (ncomp_t c=0; c<ncomp; ++c)
182 : : {
183 : : // DG(P0)
184 : 1688427 : auto mark = c*ndof;
185 : 1688427 : auto rmark = c*rdof;
186 : 1688427 : unk(e, rmark, offset) = R[mark] / L(e, mark, offset);
187 : :
188 : : // if P0P1, initialize higher dofs to 0
189 [ + + ]: 1688427 : if (rdof > ndof)
190 : : {
191 : 73085 : unk(e, rmark+1, offset) = 0.0;
192 : 73085 : unk(e, rmark+2, offset) = 0.0;
193 : 73085 : unk(e, rmark+3, offset) = 0.0;
194 : : }
195 : :
196 [ + + ]: 1688427 : if(ndof > 1) // DG(P1)
197 : : {
198 : 369127 : unk(e, rmark+1, offset) = R[mark+1] / L(e, mark+1, offset);
199 : 369127 : unk(e, rmark+2, offset) = R[mark+2] / L(e, mark+2, offset);
200 : 369127 : unk(e, rmark+3, offset) = R[mark+3] / L(e, mark+3, offset);
201 : :
202 [ + + ]: 369127 : if(ndof > 4) // DG(P2)
203 : : {
204 : 241656 : unk(e, rmark+4, offset) = R[mark+4] / L(e, mark+4, offset);
205 : 241656 : unk(e, rmark+5, offset) = R[mark+5] / L(e, mark+5, offset);
206 : 241656 : unk(e, rmark+6, offset) = R[mark+6] / L(e, mark+6, offset);
207 : 241656 : unk(e, rmark+7, offset) = R[mark+7] / L(e, mark+7, offset);
208 : 241656 : unk(e, rmark+8, offset) = R[mark+8] / L(e, mark+8, offset);
209 : 241656 : unk(e, rmark+9, offset) = R[mark+9] / L(e, mark+9, offset);
210 : : }
211 : : }
212 : : }
213 : 1076969 : }
|
