Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Box.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Initial conditions box related functionality
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Box.hpp"
14 : : #include "EOS.hpp"
15 : : #include "InciterConfig.hpp"
16 : :
17 : : namespace inciter {
18 : :
19 : : extern ctr::Config g_cfg;
20 : :
21 : : } // ::inciter
22 : :
23 : : namespace problems {
24 : :
25 : : using inciter::g_cfg;
26 : :
27 : : std::vector< std::unordered_set< std::size_t > >
28 : 2459 : boxnodes( const std::array< std::vector< tk::real >, 3 >& coord )
29 : : // *****************************************************************************
30 : : // Determine nodes that lie inside user-defined IC box(es)
31 : : //! \param[in] coord Mesh node coordinates
32 : : //! \return inbox List of nodes at which box user ICs are set for each IC box
33 : : // *****************************************************************************
34 : : {
35 : 2459 : const auto& icbox = g_cfg.get< tag::ic >();
36 : :
37 [ + + ]: 2459 : if (icbox.empty()) return {};
38 : :
39 : 3 : std::vector< std::unordered_set< std::size_t > > inbox;
40 : :
41 : 3 : const auto& x = coord[0];
42 : 3 : const auto& y = coord[1];
43 : 3 : const auto& z = coord[2];
44 : :
45 : 3 : std::size_t bcnt = 0;
46 [ + + ]: 9 : for (const auto& b : icbox) {
47 [ + - ]: 6 : inbox.emplace_back();
48 : 6 : const auto& bx = b.get< tag::x >();
49 : 6 : const auto& by = b.get< tag::y >();
50 : 6 : const auto& bz = b.get< tag::z >();
51 [ + - ]: 6 : std::vector< tk::real > box{ bx[0], bx[1], by[0], by[1], bz[0], bz[1] };
52 : :
53 : 6 : const auto eps = std::numeric_limits< tk::real >::epsilon();
54 : : // Determine which nodes lie in the IC box
55 [ + - ][ + - ]: 6 : if ( std::any_of( begin(box), end(box), [=](auto p)
56 : 6 : { return abs(p) > eps; } ) )
57 : : {
58 : 6 : std::array< tk::real, 3 > b_min{{box[0], box[2], box[4]}};
59 : 6 : std::array< tk::real, 3 > b_max{{box[1], box[3], box[5]}};
60 [ + + ]: 2952 : for (std::size_t i=0; i<x.size(); ++i) {
61 : 2946 : std::array< tk::real, 3 > node{{ x[i], y[i], z[i] }};
62 [ + + ]: 5145 : if ( node[0]>b_min[0] && node[0]<b_max[0] &&
63 [ + - ][ + - ]: 1473 : node[1]>b_min[1] && node[1]<b_max[1] &&
64 [ + + ][ + - ]: 5145 : node[2]>b_min[2] && node[2]<b_max[2] )
[ + - ][ + + ]
65 : : {
66 [ + - ]: 1473 : inbox[bcnt].insert( i );
67 : : }
68 : : }
69 : : }
70 : 6 : ++bcnt;
71 : 6 : }
72 : :
73 : 3 : return inbox;
74 : 3 : }
75 : :
76 : : void
77 : 3208319 : box( std::size_t p,
78 : : std::vector< tk::real >& u,
79 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes )
80 : : // *****************************************************************************
81 : : // Evaluate solution in user-defined IC box
82 : : //! \param[in] p Local mesh node id at which to evaluate
83 : : //! \param[in,out] u Solution to overwrite with box value
84 : : //! \param[in] boxnodes Nodes at which box user ICs are set (for each box IC)
85 : : // *****************************************************************************
86 : : {
87 : 3208319 : const auto& icbox = g_cfg.get< tag::ic >();
88 [ + + ]: 3208319 : if (icbox.empty()) return;
89 : :
90 : 1957 : auto large = std::numeric_limits< double >::max();
91 : :
92 [ + + ]: 5871 : for (std::size_t j=0; j<icbox.size(); ++j) {
93 : 3914 : const auto& b = icbox[j];
94 [ + - ][ + - ]: 3914 : if (boxnodes.size() > j && boxnodes[j].find(p) != boxnodes[j].end()) {
[ + + ][ + + ]
95 : 1957 : auto boxr = b.get< tag::ic_density >();
96 : 1957 : const auto& boxv = b.get< tag::ic_velocity >();
97 : 1957 : auto boxp = b.get< tag::ic_pressure >();
98 : 1957 : auto boxe = b.get< tag::ic_energy >();
99 : 1957 : auto boxt = b.get< tag::ic_temperature >();
100 : :
101 : 1957 : tk::real r = 0.0, ru = 0.0, rv = 0.0, rw = 0.0, re = 0.0;
102 [ + - ][ + - ]: 1957 : if (std::abs(boxr - large) > 1.0e-12 && boxr > 0.0) {
[ + - ]
103 : 1957 : r = boxr;
104 : : }
105 [ + - ]: 3914 : if (std::abs(boxv[0] - large) > 1.0e-12 &&
106 [ + - ][ + - ]: 3914 : std::abs(boxv[1] - large) > 1.0e-12 &&
[ + - ]
107 : 1957 : std::abs(boxv[2] - large) > 1.0e-12)
108 : : {
109 : 1957 : ru = r * boxv[0];
110 : 1957 : rv = r * boxv[1];
111 : 1957 : rw = r * boxv[2];
112 : : }
113 [ + - ][ + - ]: 1957 : if (std::abs(boxp - large) > 1.0e-12 && boxp> 0.0) {
[ + - ]
114 : 1957 : re = eos::totalenergy( r, ru/r, rv/r, rw/r, boxp);
115 : : }
116 [ - + ][ - - ]: 1957 : if (std::abs(boxe - large) > 1.0e-12 && boxe > 0.0) {
[ - + ]
117 : 0 : auto ux = ru/r, uy = rv/r, uz = rw/r;
118 : 0 : auto ke = 0.5*(ux*ux + uy*uy + uz*uz);
119 : 0 : re = r * (boxe + ke);
120 : : }
121 [ - + ][ - - ]: 1957 : if (std::abs(boxt - large) > 1.0e-12 && boxt > 0.0) {
[ - + ]
122 : 0 : auto cv = g_cfg.get< tag::mat_spec_heat_const_vol >();
123 [ - - ][ - - ]: 0 : if (std::abs(cv - large) > 1.0e-12 && cv > 0.0) {
[ - - ]
124 : 0 : re = r * boxt * cv;
125 : : }
126 : : }
127 : 1957 : u[0] = r;
128 : 1957 : u[1] = ru;
129 : 1957 : u[2] = rv;
130 : 1957 : u[3] = rw;
131 : 1957 : u[4] = re;
132 : : }
133 : : }
134 : : }
135 : :
136 : : } // problems::
|
