Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/CompFlow/Physics/CGNavierStokes.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 Physics policy for the Navier-Stokes equation using continuous
9 : : Galerkin discretization
10 : : \details This file defines a Physics policy class for solving the
11 : : compressible single-material viscous flow equations using continuous
12 : : Galerkin discretization, defined in PDE/CompFlow/CGCompFlow.h. The class
13 : : defined here is used to configure the behavior of CGCompFlow. See
14 : : PDE/CompFlow/Physics/CG.h for general requirements on Physics policy classes
15 : : for CGCompFlow.
16 : : */
17 : : // *****************************************************************************
18 : :
19 : : #include "CGNavierStokes.hpp"
20 : : #include "EoS/EoS.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : :
26 : : } // inciter::
27 : :
28 : : using inciter::cg::CompFlowPhysicsNavierStokes;
29 : :
30 : : void
31 : 0 : CompFlowPhysicsNavierStokes::viscousRhs(
32 : : tk::real dt,
33 : : tk::real J,
34 : : const std::array< std::size_t, 4 >& N,
35 : : const std::array< std::array< tk::real, 3 >, 4 >& grad,
36 : : const std::array< std::array< tk::real, 4 >, 5 >& u,
37 : : const std::array< const tk::real*, 5 >& r,
38 : : tk::Fields& R ) const
39 : : // *****************************************************************************
40 : : // Add viscous stress contribution to momentum and energy rhs
41 : : //! \param[in] dt Size of time step
42 : : //! \param[in] J Element Jacobi determinant
43 : : //! \param[in] N Element node indices
44 : : //! \param[in] grad Shape function derivatives, nnode*ndim [4][3]
45 : : //! \param[in] u Solution at element nodes at recent time step
46 : : //! \param[in] r Pointers to right hand side at component and offset
47 : : //! \param[in,out] R Right-hand side vector contributing to
48 : : // *****************************************************************************
49 : : {
50 : : // dynamic viscosity
51 : : auto mu_d = mu< tag::compflow >(0);
52 : :
53 : : // add deviatoric viscous stress contribution to momentum rhs
54 : 0 : auto c = dt * J/6.0 * mu_d;
55 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
56 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
57 [ - - ]: 0 : for (std::size_t k=0; k<4; ++k) {
58 : 0 : R.var(r[i+1],N[0]) -= c * grad[0][j]*(grad[k][j]*u[i+1][k] +
59 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
60 : 0 : R.var(r[i+1],N[1]) -= c * grad[1][j]*(grad[k][j]*u[i+1][k] +
61 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
62 : 0 : R.var(r[i+1],N[2]) -= c * grad[2][j]*(grad[k][j]*u[i+1][k] +
63 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
64 : 0 : R.var(r[i+1],N[3]) -= c * grad[3][j]*(grad[k][j]*u[i+1][k] +
65 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
66 : : }
67 : : // add isotropic viscous stress contribution to momentum rhs
68 : 0 : c = dt * J/6.0 * mu_d * 2.0/3.0;
69 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
70 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
71 [ - - ]: 0 : for (std::size_t k=0; k<4; ++k) {
72 : 0 : R.var(r[i+1],N[0]) += c * grad[0][i]*grad[k][j]*u[j+1][k]/u[0][k];
73 : 0 : R.var(r[i+1],N[1]) += c * grad[1][i]*grad[k][j]*u[j+1][k]/u[0][k];
74 : 0 : R.var(r[i+1],N[2]) += c * grad[2][i]*grad[k][j]*u[j+1][k]/u[0][k];
75 : 0 : R.var(r[i+1],N[3]) += c * grad[3][i]*grad[k][j]*u[j+1][k]/u[0][k];
76 : : }
77 : : // add deviatoric viscous stress contribution to energy rhs
78 : 0 : c = dt * J/24.0 * mu_d;
79 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
80 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
81 [ - - ]: 0 : for (std::size_t k=0; k<4; ++k) {
82 : 0 : R.var(r[4],N[0]) -= c * u[i+1][k]/u[0][k] *
83 : 0 : grad[0][j]*(grad[k][j]*u[i+1][k] +
84 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
85 : 0 : R.var(r[4],N[1]) -= c * u[i+1][k]/u[0][k] *
86 : 0 : grad[1][j]*(grad[k][j]*u[i+1][k] +
87 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
88 : 0 : R.var(r[4],N[2]) -= c * u[i+1][k]/u[0][k] *
89 : 0 : grad[2][j]*(grad[k][j]*u[i+1][k] +
90 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
91 : 0 : R.var(r[4],N[3]) -= c * u[i+1][k]/u[0][k] *
92 : 0 : grad[3][j]*(grad[k][j]*u[i+1][k] +
93 : 0 : grad[k][i]*u[j+1][k])/u[0][k];
94 : : }
95 : : // add isotropic viscous stress contribution to energy rhs
96 : 0 : c = dt * J/24.0 * mu_d * 2.0/3.0;
97 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
98 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
99 [ - - ]: 0 : for (std::size_t k=0; k<4; ++k) {
100 : 0 : R.var(r[4],N[0]) += c * u[i+1][k]/u[0][k] *
101 : 0 : grad[0][i]*grad[k][j]*u[j+1][k]/u[0][k];
102 : 0 : R.var(r[4],N[1]) += c * u[i+1][k]/u[0][k] *
103 : 0 : grad[1][i]*grad[k][j]*u[j+1][k]/u[0][k];
104 : 0 : R.var(r[4],N[2]) += c * u[i+1][k]/u[0][k] *
105 : 0 : grad[2][i]*grad[k][j]*u[j+1][k]/u[0][k];
106 : 0 : R.var(r[4],N[3]) += c * u[i+1][k]/u[0][k] *
107 : 0 : grad[3][i]*grad[k][j]*u[j+1][k]/u[0][k];
108 : : }
109 : 0 : }
110 : :
111 : : tk::real
112 : 0 : CompFlowPhysicsNavierStokes::viscous_dt(
113 : : tk::real L,
114 : : const std::array< std::array< tk::real, 4 >, 5 >& u ) const
115 : : // *****************************************************************************
116 : : // Compute the minimum time step size based on the viscous force
117 : : //! \param[in] L Characteristic length scale
118 : : //! \param[in] u Solution at element nodes at recent time step
119 : : //! \return Minimum time step size based on viscous force
120 : : // *****************************************************************************
121 : : {
122 : : // dynamic viscosity
123 : : auto mu_d = mu< tag::compflow >(0);
124 : :
125 : : // compute the minimum viscous time step size across the four nodes
126 : : tk::real mindt = std::numeric_limits< tk::real >::max();
127 [ - - ]: 0 : for (std::size_t j=0; j<4; ++j) {
128 : : auto& r = u[0][j]; // rho
129 : 0 : auto dt = L * L * r / (2.0*mu_d); // dt ~ dx^2/nu = dx^2*rho/(2mu)
130 [ - - ]: 0 : if (dt < mindt) mindt = dt;
131 : : }
132 : 0 : return mindt;
133 : : }
134 : :
135 : : void
136 : 0 : CompFlowPhysicsNavierStokes::conductRhs(
137 : : tk::real dt,
138 : : tk::real J,
139 : : const std::array< std::size_t, 4 >& N,
140 : : const std::array< std::array< tk::real, 3 >, 4 >& grad,
141 : : const std::array< std::array< tk::real, 4 >, 5 >& u,
142 : : const std::array< const tk::real*, 5 >& r,
143 : : tk::Fields& R ) const
144 : : // *****************************************************************************
145 : : //! Add heat conduction contribution to the energy rhs
146 : : //! \param[in] dt Size of time step
147 : : //! \param[in] J Element Jacobi determinant
148 : : //! \param[in] N Element node indices
149 : : //! \param[in] grad Shape function derivatives, nnode*ndim [4][3]
150 : : //! \param[in] u Solution at element nodes at recent time step
151 : : //! \param[in] r Pointers to right hand side at component and offset
152 : : //! \param[in,out] R Right-hand side vector contributing to
153 : : // *****************************************************************************
154 : : {
155 : : // specific heat at constant volume
156 : : auto c_v = cv< tag::compflow >(0);
157 : : // thermal conductivity
158 : : auto kc = k< tag::compflow >(0);
159 : :
160 : : // compute temperature
161 : : std::array< tk::real, 4 > T;
162 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i)
163 : 0 : T[i] = c_v*(u[4][i] - (u[1][i]*u[1][i] +
164 : 0 : u[2][i]*u[2][i] +
165 : 0 : u[3][i]*u[3][i])/2.0/u[0][i]) / u[0][i];
166 : : // add heat conduction contribution to energy rhs
167 : 0 : auto c = dt * J/24.0 * kc;
168 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
169 [ - - ]: 0 : for (std::size_t k=0; k<4; ++k) {
170 : 0 : R.var(r[4],N[0]) += c * grad[k][i] * T[k];
171 : 0 : R.var(r[4],N[1]) += c * grad[k][i] * T[k];
172 : 0 : R.var(r[4],N[2]) += c * grad[k][i] * T[k];
173 : 0 : R.var(r[4],N[3]) += c * grad[k][i] * T[k];
174 : : }
175 : 0 : }
176 : :
177 : : tk::real
178 : 0 : CompFlowPhysicsNavierStokes::conduct_dt(
179 : : tk::real L,
180 : : tk::real g,
181 : : const std::array< std::array< tk::real, 4 >, 5 >& u ) const
182 : : // *****************************************************************************
183 : : //! Compute the minimum time step size based on thermal diffusion
184 : : //! \param[in] L Characteristic length scale
185 : : //! \param[in] g Ratio of specific heats
186 : : //! \param[in] u Solution at element nodes at recent time step
187 : : //! \return Minimum time step size based on thermal diffusion
188 : : // *****************************************************************************
189 : : {
190 : : // specific heat at constant volume
191 : : auto c_v = cv< tag::compflow >(0);
192 : : // thermal conductivity
193 : : auto kc = k< tag::compflow >(0);
194 : : // specific heat at constant pressure
195 : 0 : auto cp = g * c_v;
196 : :
197 : : // compute the minimum conduction time step size across the four nodes
198 : : tk::real mindt = std::numeric_limits< tk::real >::max();
199 [ - - ]: 0 : for (std::size_t j=0; j<4; ++j) {
200 : : auto& r = u[0][j]; // rho
201 : 0 : auto dt = L * L * r * cp / kc; // dt ~ dx^2/alpha = dx^2*rho*C_p/kc
202 [ - - ]: 0 : if (dt < mindt) mindt = dt;
203 : : }
204 : 0 : return mindt;
205 : : }
|
