Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Kozak.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 KozCG: Taylor-Galerkin, FCT, element-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Vector.hpp"
14 : : #include "EOS.hpp"
15 : : #include "Kozak.hpp"
16 : : #include "Problems.hpp"
17 : : #include "InciterConfig.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : : extern ctr::Config g_cfg;
22 : :
23 : : } // ::inciter
24 : :
25 : : namespace kozak {
26 : :
27 : : using inciter::g_cfg;
28 : :
29 : : void
30 : 9006 : rhs( const std::vector< std::size_t >& inpoel,
31 : : const std::array< std::vector< tk::real >, 3 >& coord,
32 : : tk::real t,
33 : : tk::real dt,
34 : : const std::vector< tk::real >& tp,
35 : : const std::vector< tk::real >& dtp,
36 : : const tk::Fields& U,
37 : : tk::Fields& R )
38 : : // *****************************************************************************
39 : : // Compute right hand side
40 : : //! \param[in] inpoel Tetrahedron connectivity
41 : : //! \param[in] coord Mesh node coordinates
42 : : //! \param[in] U Unknowns/solution vector in mesh nodes
43 : : //! \param[in] t Physical time
44 : : //! \param[in] dt Physical time size
45 : : //! \param[in] tp Phisical time step size for each mesh node (if steady state)
46 : : //! \param[in] dtp Time step size for each mesh node (if steady state)
47 : : //! \param[in,out] R Right-hand side vector computed
48 : : // *****************************************************************************
49 : : {
50 [ - + ][ - - ]: 9006 : Assert( U.nunk() == coord[0].size(), "Size mismatch" );
[ - - ][ - - ]
51 [ - + ][ - - ]: 9006 : Assert( R.nunk() == coord[0].size(), "Size mismatch" );
[ - - ][ - - ]
52 : :
53 : : // zero right hand side for all components
54 [ + - ]: 9006 : R.fill( 0.0 );
55 : :
56 : 9006 : const auto steady = g_cfg.get< tag::steady >();
57 : 9006 : const auto ncomp = U.nprop();
58 : 9006 : const auto& x = coord[0];
59 : 9006 : const auto& y = coord[1];
60 : 9006 : const auto& z = coord[2];
61 : :
62 [ + - ]: 9006 : auto src = problems::SRC();
63 : :
64 : : #if defined(__clang__)
65 : : #pragma clang diagnostic push
66 : : #pragma clang diagnostic ignored "-Wvla"
67 : : #pragma clang diagnostic ignored "-Wvla-extension"
68 : : #elif defined(STRICT_GNUC)
69 : : #pragma GCC diagnostic push
70 : : #pragma GCC diagnostic ignored "-Wvla"
71 : : #endif
72 : :
73 [ + + ]: 932306 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
74 : :
75 : : // Element gradients and Jacobian
76 : 1846600 : const auto N = inpoel.data() + e*4;
77 : : const std::array< tk::real, 3 >
78 : 923300 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
79 : 923300 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
80 : 923300 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
81 : 923300 : const auto J = tk::triple( ba, ca, da );
82 : : std::array< std::array< tk::real, 3 >, 4 > grad;
83 : 923300 : grad[1] = tk::cross( ca, da );
84 : 923300 : grad[2] = tk::cross( da, ba );
85 : 923300 : grad[3] = tk::cross( ba, ca );
86 [ + + ]: 3693200 : for (std::size_t i=0; i<3; ++i) {
87 : 2769900 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
88 : : }
89 : :
90 : : // Taylor-Galerkin first half step
91 : :
92 : : tk::real p[4];
93 [ + + ]: 4616500 : for (std::size_t a=0; a<4; ++a) {
94 [ + - ]: 3693200 : auto r = U(N[a],0);
95 [ + - ]: 3693200 : auto ru = U(N[a],1);
96 [ + - ]: 3693200 : auto rv = U(N[a],2);
97 [ + - ]: 3693200 : auto rw = U(N[a],3);
98 [ + - ]: 3693200 : p[a] = eos::pressure( U(N[a],4) - 0.5*(ru*ru + rv*rv + rw*rw)/r );
99 : : }
100 : :
101 : 923300 : tk::real ue[ncomp];
102 [ + + ]: 5758380 : for (std::size_t c=0; c<ncomp; ++c) {
103 [ + - ][ + - ]: 4835080 : ue[c] = (U(N[0],c) + U(N[1],c) + U(N[2],c) + U(N[3],c))/4.0;
[ + - ][ + - ]
104 : : }
105 : :
106 [ - + ]: 923300 : if (steady) dt = (dtp[N[0]] + dtp[N[1]] + dtp[N[2]] + dtp[N[3]])/4.0;
107 : :
108 : 923300 : auto coef = dt/J/2.0;
109 [ + + ]: 3693200 : for (std::size_t j=0; j<3; ++j) {
110 [ + + ]: 13849500 : for (std::size_t a=0; a<4; ++a) {
111 : 11079600 : auto cg = coef * grad[a][j];
112 [ + - ][ + - ]: 11079600 : auto uj = U(N[a],j+1) / U(N[a],0);
113 [ + - ]: 11079600 : ue[0] -= cg * U(N[a],j+1);
114 [ + - ]: 11079600 : ue[1] -= cg * U(N[a],1) * uj;
115 [ + - ]: 11079600 : ue[2] -= cg * U(N[a],2) * uj;
116 [ + - ]: 11079600 : ue[3] -= cg * U(N[a],3) * uj;
117 : 11079600 : ue[j+1] -= cg * p[a];
118 [ + - ]: 11079600 : ue[4] -= cg * (U(N[a],4) + p[a]) * uj;
119 [ + + ]: 13702560 : for (std::size_t c=5; c<ncomp; ++c) {
120 [ + - ]: 2622960 : ue[c] -= cg * U(N[a],c) * uj;
121 : : }
122 : : }
123 : : }
124 : :
125 [ + + ]: 923300 : if (src) {
126 : 685780 : coef = dt/8.0;
127 [ - + ]: 685780 : if (steady) t = (tp[N[0]] + tp[N[1]] + tp[N[2]] + tp[N[3]])/4.0;
128 [ + + ]: 3428900 : for (std::size_t a=0; a<4; ++a) {
129 [ + - ]: 2743120 : auto s = src( x[N[a]], y[N[a]], z[N[a]], t );
130 [ + + ]: 17333040 : for (std::size_t c=0; c<ncomp; ++c) {
131 : 14589920 : ue[c] += coef * s[c];
132 : : }
133 : 2743120 : }
134 : : }
135 : :
136 : : // Taylor-Galerkin: second half step
137 : :
138 : 923300 : auto r = ue[0];
139 : 923300 : auto ru = ue[1];
140 : 923300 : auto rv = ue[2];
141 : 923300 : auto rw = ue[3];
142 : 923300 : auto pr = eos::pressure( ue[4] - 0.5*(ru*ru + rv*rv + rw*rw)/r );
143 : :
144 : 923300 : coef = 1.0/6.0;
145 [ + + ]: 3693200 : for (std::size_t j=0; j<3; ++j) {
146 : 2769900 : auto uj = ue[j+1] / ue[0];
147 [ + + ]: 13849500 : for (std::size_t a=0; a<4; ++a) {
148 : 11079600 : auto cg = coef * grad[a][j];
149 [ + - ]: 11079600 : R(N[a],0) += cg * ue[j+1];
150 [ + - ]: 11079600 : R(N[a],1) += cg * ue[1] * uj;
151 [ + - ]: 11079600 : R(N[a],2) += cg * ue[2] * uj;
152 [ + - ]: 11079600 : R(N[a],3) += cg * ue[3] * uj;
153 [ + - ]: 11079600 : R(N[a],j+1) += cg * pr;
154 [ + - ]: 11079600 : R(N[a],4) += cg * (ue[4] + pr) * uj;
155 [ + + ]: 13702560 : for (std::size_t c=5; c<ncomp; ++c) {
156 [ + - ]: 2622960 : R(N[a],c) += cg * ue[c] * uj;
157 : : }
158 : : }
159 : : }
160 : :
161 [ + + ]: 923300 : if (src) {
162 : 685780 : auto xe = (x[N[0]] + x[N[1]] + x[N[2]] + x[N[3]])/4.0;
163 : 685780 : auto ye = (y[N[0]] + y[N[1]] + y[N[2]] + y[N[3]])/4.0;
164 : 685780 : auto ze = (z[N[0]] + z[N[1]] + z[N[2]] + z[N[3]])/4.0;
165 [ + - ]: 685780 : auto se = src( xe, ye, ze, t+dt/2.0 );
166 : 685780 : coef = J/24.0;
167 [ + + ]: 3428900 : for (std::size_t a=0; a<4; ++a) {
168 [ + + ]: 17333040 : for (std::size_t c=0; c<ncomp; ++c) {
169 [ + - ]: 14589920 : R(N[a],c) += coef * se[c];
170 : : }
171 : : }
172 : 685780 : }
173 : 923300 : }
174 : :
175 : : #if defined(__clang__)
176 : : #pragma clang diagnostic pop
177 : : #elif defined(STRICT_GNUC)
178 : : #pragma GCC diagnostic pop
179 : : #endif
180 : 9006 : }
181 : :
182 : : } // kozak::
|
