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