Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Initialize.hpp
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 : : #ifndef Initialize_h
15 : : #define Initialize_h
16 : :
17 : : #include "Basis.hpp"
18 : : #include "Fields.hpp"
19 : : #include "UnsMesh.hpp"
20 : : #include "FunctionPrototypes.hpp"
21 : : #include "Inciter/InputDeck/InputDeck.hpp"
22 : : #include "EoS/EOS.hpp"
23 : :
24 : : namespace inciter {
25 : :
26 : : extern ctr::InputDeck g_inputdeck;
27 : :
28 : : } // inciter::
29 : :
30 : : namespace tk {
31 : :
32 : : //! Initalize a PDE system for DG by projecting the exact solution
33 : : //! in the DG solution space
34 : : void
35 : : initialize( ncomp_t ncomp,
36 : : const std::vector< inciter::EOS >& mat_blk,
37 : : const Fields& L,
38 : : const std::vector< std::size_t >& inpoel,
39 : : const UnsMesh::Coords& coord,
40 : : const InitializeFn& solution,
41 : : Fields& unk,
42 : : real t,
43 : : const std::size_t nielem );
44 : :
45 : : //! Update the rhs by adding the initial analytical solution term
46 : : void
47 : : update_rhs( ncomp_t ncomp,
48 : : const std::size_t ndof,
49 : : const tk::real wt,
50 : : const std::vector< tk::real >& B,
51 : : const std::vector< tk::real >& s,
52 : : std::vector< tk::real >& R );
53 : :
54 : :
55 : : //! Compute the initial conditions
56 : : void
57 : : eval_init( ncomp_t ncomp,
58 : : const std::size_t ndof,
59 : : const std::size_t rdof,
60 : : const std::size_t e,
61 : : const std::vector< tk::real >& R,
62 : : const Fields& L,
63 : : Fields& unk );
64 : :
65 : : template< class eq >
66 : : void
67 : 575 : BoxElems(
68 : : const tk::Fields& geoElem,
69 : : std::size_t nielem,
70 : : std::vector< std::unordered_set< std::size_t > >& inbox )
71 : : // *****************************************************************************
72 : : //! Determine if elements lie inside user defined IC boxes
73 : : //! \tparam eq Equation type to operate on, e.g., tag::compflow, tag::multimat
74 : : //! \param[in] geoElem Element geometry array
75 : : //! \param[in] nielem Number of internal elements
76 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
77 : : //! each IC box
78 : : // *****************************************************************************
79 : : {
80 : : // Detect if user has configured IC boxes
81 : 575 : const auto& icbox = inciter::g_inputdeck.get<tag::ic, tag::box>();
82 [ + + ]: 575 : if (!icbox.empty()) {
83 : 5 : std::size_t bcnt = 0;
84 [ + + ]: 10 : for (const auto& b : icbox) { // for all boxes for this eq
85 [ + - ]: 5 : inbox.emplace_back();
86 [ + - ]: 5 : std::vector< tk::real > box
87 : 5 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
88 : 5 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
89 : 5 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
90 : :
91 : : // Determine orientation of box
92 : 15 : std::array< tk::real, 3 > b_orientn{{
93 : 5 : b.template get< tag::orientation >()[0],
94 : 5 : b.template get< tag::orientation >()[1],
95 : 5 : b.template get< tag::orientation >()[2] }};
96 : 5 : std::array< tk::real, 3 > b_centroid{{ 0.5*(box[0]+box[1]),
97 : 5 : 0.5*(box[2]+box[3]), 0.5*(box[4]+box[5]) }};
98 : :
99 : 5 : const auto eps = std::numeric_limits< tk::real >::epsilon();
100 : : // Determine which elements lie in the IC box
101 [ + - ][ + - ]: 5 : if ( std::any_of( begin(box), end(box),
102 : 5 : [=](auto p){ return abs(p) > eps; } ) )
103 : : {
104 : : // Transform box to reference space
105 : 5 : std::array< tk::real, 3 > b_min{{box[0], box[2], box[4]}};
106 : 5 : std::array< tk::real, 3 > b_max{{box[1], box[3], box[5]}};
107 : 5 : tk::movePoint(b_centroid, b_min);
108 : 5 : tk::movePoint(b_centroid, b_max);
109 : :
110 [ + + ]: 2417 : for (ncomp_t e=0; e<nielem; ++e) {
111 [ + - ]: 2412 : auto x = geoElem(e,1);
112 [ + - ]: 2412 : auto y = geoElem(e,2);
113 [ + - ]: 2412 : auto z = geoElem(e,3);
114 : 2412 : std::array< tk::real, 3 > node{{ x, y, z }};
115 : : // Transform node to reference space of box
116 : 2412 : tk::movePoint(b_centroid, node);
117 : 2412 : tk::rotatePoint({{-b_orientn[0], -b_orientn[1], -b_orientn[2]}},
118 : : node);
119 [ + + ]: 4824 : if ( node[0]>b_min[0] && node[0]<b_max[0] &&
120 [ + - ][ + - ]: 1046 : node[1]>b_min[1] && node[1]<b_max[1] &&
121 [ + - ][ + - ]: 4824 : node[2]>b_min[2] && node[2]<b_max[2] )
[ + - ][ + + ]
122 : : {
123 [ + - ]: 1046 : inbox[bcnt].insert( e );
124 : : }
125 : : }
126 : : }
127 : 5 : ++bcnt;
128 : : }
129 : : }
130 : 575 : }
131 : :
132 : : } // tk::
133 : :
134 : : #endif // Initialize_h
|
