Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PCNSFVKT.h"
11 : #include "NS.h"
12 : #include "SinglePhaseFluidProperties.h"
13 : #include "Limiter.h"
14 :
15 : using namespace Moose::FV;
16 :
17 : registerMooseObject("NavierStokesApp", PCNSFVKT);
18 :
19 : InputParameters
20 1535 : PCNSFVKT::validParams()
21 : {
22 1535 : InputParameters params = FVFluxKernel::validParams();
23 1535 : params.addClassDescription("Computes the residual of advective term using finite volume method.");
24 1535 : params.addRequiredParam<UserObjectName>(NS::fluid, "Fluid userobject");
25 3070 : MooseEnum eqn("mass momentum energy scalar");
26 3070 : params.addRequiredParam<MooseEnum>("eqn", eqn, "The equation you're solving.");
27 3070 : MooseEnum momentum_component("x=0 y=1 z=2");
28 3070 : params.addParam<MooseEnum>("momentum_component",
29 : momentum_component,
30 : "The component of the momentum equation that this kernel applies to.");
31 3070 : params.addParam<MooseEnum>(
32 : "limiter", moose_limiter_type, "The limiter to apply during interpolation.");
33 1535 : params.set<unsigned short>("ghost_layers") = 2;
34 3070 : params.addParam<bool>(
35 : "knp_for_omega",
36 3070 : true,
37 : "Whether to use the Kurganov, Noelle, and Petrova method to compute the omega parameter for "
38 : "stabilization. If false, then the Kurganov-Tadmor method will be used.");
39 1535 : return params;
40 1535 : }
41 :
42 779 : PCNSFVKT::PCNSFVKT(const InputParameters & params)
43 : : FVFluxKernel(params),
44 1558 : _fluid(getUserObject<SinglePhaseFluidProperties>(NS::fluid)),
45 :
46 779 : _sup_vel_x_elem(getADMaterialProperty<Real>(NS::superficial_velocity_x)),
47 779 : _sup_vel_x_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_x)),
48 779 : _grad_sup_vel_x_elem(
49 779 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_x))),
50 779 : _grad_sup_vel_x_neighbor(
51 779 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_x))),
52 779 : _sup_vel_y_elem(getADMaterialProperty<Real>(NS::superficial_velocity_y)),
53 779 : _sup_vel_y_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_y)),
54 779 : _grad_sup_vel_y_elem(
55 779 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_y))),
56 779 : _grad_sup_vel_y_neighbor(
57 779 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_y))),
58 779 : _sup_vel_z_elem(getADMaterialProperty<Real>(NS::superficial_velocity_z)),
59 779 : _sup_vel_z_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_z)),
60 779 : _grad_sup_vel_z_elem(
61 779 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_z))),
62 779 : _grad_sup_vel_z_neighbor(
63 779 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_z))),
64 :
65 779 : _T_fluid_elem(getADMaterialProperty<Real>(NS::T_fluid)),
66 779 : _T_fluid_neighbor(getNeighborADMaterialProperty<Real>(NS::T_fluid)),
67 779 : _grad_T_fluid_elem(getADMaterialProperty<RealVectorValue>(NS::grad(NS::T_fluid))),
68 779 : _grad_T_fluid_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::T_fluid))),
69 779 : _pressure_elem(getADMaterialProperty<Real>(NS::pressure)),
70 779 : _pressure_neighbor(getNeighborADMaterialProperty<Real>(NS::pressure)),
71 779 : _grad_pressure_elem(getADMaterialProperty<RealVectorValue>(NS::grad(NS::pressure))),
72 779 : _grad_pressure_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::pressure))),
73 779 : _eps_elem(getMaterialProperty<Real>(NS::porosity)),
74 779 : _eps_neighbor(getNeighborMaterialProperty<Real>(NS::porosity)),
75 1558 : _eqn(getParam<MooseEnum>("eqn")),
76 1558 : _index(getParam<MooseEnum>("momentum_component")),
77 779 : _scalar_elem(_u_elem),
78 779 : _scalar_neighbor(_u_neighbor),
79 779 : _grad_scalar_elem((_eqn == "scalar") ? &_var.adGradSln() : nullptr),
80 779 : _grad_scalar_neighbor((_eqn == "scalar") ? &_var.adGradSlnNeighbor() : nullptr),
81 1558 : _limiter(Limiter<ADReal>::build(LimiterType(int(getParam<MooseEnum>("limiter"))))),
82 2337 : _knp_for_omega(getParam<bool>("knp_for_omega"))
83 : {
84 1619 : if ((_eqn == "momentum") && !isParamValid("momentum_component"))
85 0 : paramError("eqn",
86 : "If 'momentum' is specified for 'eqn', then you must provide a parameter "
87 : "value for 'momentum_component'");
88 779 : }
89 :
90 : std::pair<ADReal, ADReal>
91 8839738 : PCNSFVKT::computeAlphaAndOmega(const ADReal & u_elem_normal,
92 : const ADReal & u_neighbor_normal,
93 : const ADReal & c_elem,
94 : const ADReal & c_neighbor) const
95 : {
96 : // Equations from Greenshields sans multiplication by area (which will be done in
97 : // computeResidual/Jacobian
98 : const auto psi_elem =
99 8839738 : std::max({c_elem + u_elem_normal, c_neighbor + u_neighbor_normal, ADReal(0)});
100 : const auto psi_neighbor =
101 17679476 : std::max({c_elem - u_elem_normal, c_neighbor - u_neighbor_normal, ADReal(0)});
102 17603300 : auto alpha = _knp_for_omega ? psi_elem / (psi_elem + psi_neighbor) : ADReal(0.5);
103 52733724 : auto omega = _knp_for_omega ? alpha * (1 - alpha) * (psi_elem + psi_neighbor)
104 : : alpha * std::max(psi_elem, psi_neighbor);
105 :
106 : // Do this to avoid new nonzero mallocs
107 17679476 : const auto dummy_quant = 0 * (c_elem + u_elem_normal + c_neighbor + u_neighbor_normal);
108 :
109 8839738 : alpha += dummy_quant;
110 8839738 : omega += dummy_quant;
111 8839738 : return std::make_pair(std::move(alpha), std::move(omega));
112 17679476 : }
113 :
114 : ADReal
115 8839738 : PCNSFVKT::computeFaceFlux(const ADReal & alpha,
116 : const ADReal & omega,
117 : const ADReal & sup_vel_elem_normal,
118 : const ADReal & sup_vel_neighbor_normal,
119 : const ADReal & adv_quant_elem,
120 : const ADReal & adv_quant_neighbor)
121 : {
122 8839738 : return alpha * (sup_vel_elem_normal * adv_quant_elem) +
123 8839738 : (1 - alpha) * sup_vel_neighbor_normal * adv_quant_neighbor -
124 8839738 : omega * (adv_quant_neighbor - adv_quant_elem);
125 : }
126 :
127 : ADReal
128 8839738 : PCNSFVKT::computeQpResidual()
129 : {
130 : // Perform primitive interpolations
131 : const auto pressure_elem = interpolate(*_limiter,
132 8839738 : _pressure_elem[_qp],
133 8839738 : _pressure_neighbor[_qp],
134 8839738 : &_grad_pressure_elem[_qp],
135 8839738 : *_face_info,
136 8839738 : /*elem_is_up=*/true);
137 : const auto pressure_neighbor = interpolate(*_limiter,
138 8839738 : _pressure_neighbor[_qp],
139 8839738 : _pressure_elem[_qp],
140 8839738 : &_grad_pressure_neighbor[_qp],
141 8839738 : *_face_info,
142 8839738 : /*elem_is_up=*/false);
143 : const auto T_fluid_elem = interpolate(*_limiter,
144 8839738 : _T_fluid_elem[_qp],
145 8839738 : _T_fluid_neighbor[_qp],
146 8839738 : &_grad_T_fluid_elem[_qp],
147 8839738 : *_face_info,
148 8839738 : /*elem_is_up=*/true);
149 : const auto T_fluid_neighbor = interpolate(*_limiter,
150 8839738 : _T_fluid_neighbor[_qp],
151 8839738 : _T_fluid_elem[_qp],
152 8839738 : &_grad_T_fluid_neighbor[_qp],
153 8839738 : *_face_info,
154 8839738 : /*elem_is_up=*/false);
155 : const auto sup_vel_x_elem = interpolate(*_limiter,
156 8839738 : _sup_vel_x_elem[_qp],
157 8839738 : _sup_vel_x_neighbor[_qp],
158 8839738 : &_grad_sup_vel_x_elem[_qp],
159 8839738 : *_face_info,
160 8839738 : /*elem_is_up=*/true);
161 : const auto sup_vel_x_neighbor = interpolate(*_limiter,
162 8839738 : _sup_vel_x_neighbor[_qp],
163 8839738 : _sup_vel_x_elem[_qp],
164 8839738 : &_grad_sup_vel_x_neighbor[_qp],
165 8839738 : *_face_info,
166 8839738 : /*elem_is_up=*/false);
167 : const auto sup_vel_y_elem = interpolate(*_limiter,
168 8839738 : _sup_vel_y_elem[_qp],
169 8839738 : _sup_vel_y_neighbor[_qp],
170 8839738 : &_grad_sup_vel_y_elem[_qp],
171 8839738 : *_face_info,
172 8839738 : /*elem_is_up=*/true);
173 : const auto sup_vel_y_neighbor = interpolate(*_limiter,
174 8839738 : _sup_vel_y_neighbor[_qp],
175 8839738 : _sup_vel_y_elem[_qp],
176 8839738 : &_grad_sup_vel_y_neighbor[_qp],
177 8839738 : *_face_info,
178 8839738 : /*elem_is_up=*/false);
179 : const auto sup_vel_z_elem = interpolate(*_limiter,
180 8839738 : _sup_vel_z_elem[_qp],
181 8839738 : _sup_vel_z_neighbor[_qp],
182 8839738 : &_grad_sup_vel_z_elem[_qp],
183 8839738 : *_face_info,
184 8839738 : /*elem_is_up=*/true);
185 : const auto sup_vel_z_neighbor = interpolate(*_limiter,
186 8839738 : _sup_vel_z_neighbor[_qp],
187 8839738 : _sup_vel_z_elem[_qp],
188 8839738 : &_grad_sup_vel_z_neighbor[_qp],
189 8839738 : *_face_info,
190 8839738 : /*elem_is_up=*/false);
191 :
192 : const auto sup_vel_elem = VectorValue<ADReal>(sup_vel_x_elem, sup_vel_y_elem, sup_vel_z_elem);
193 8839738 : const auto u_elem = sup_vel_elem / _eps_elem[_qp];
194 8839738 : const auto rho_elem = _fluid.rho_from_p_T(pressure_elem, T_fluid_elem);
195 8839738 : const auto specific_volume_elem = 1. / rho_elem;
196 8839738 : const auto e_elem = _fluid.e_from_T_v(T_fluid_elem, specific_volume_elem);
197 : const auto sup_vel_neighbor =
198 : VectorValue<ADReal>(sup_vel_x_neighbor, sup_vel_y_neighbor, sup_vel_z_neighbor);
199 8839738 : const auto u_neighbor = sup_vel_neighbor / _eps_neighbor[_qp];
200 8839738 : const auto rho_neighbor = _fluid.rho_from_p_T(pressure_neighbor, T_fluid_neighbor);
201 8839738 : const auto specific_volume_neighbor = 1. / rho_neighbor;
202 8839738 : const auto e_neighbor = _fluid.e_from_T_v(T_fluid_neighbor, specific_volume_neighbor);
203 :
204 8839738 : const auto c_elem = _fluid.c_from_v_e(specific_volume_elem, e_elem);
205 8839738 : const auto c_neighbor = _fluid.c_from_v_e(specific_volume_neighbor, e_neighbor);
206 :
207 8839738 : const auto sup_vel_elem_normal = sup_vel_elem * _face_info->normal();
208 8839738 : const auto sup_vel_neighbor_normal = sup_vel_neighbor * _face_info->normal();
209 8839738 : const auto u_elem_normal = u_elem * _face_info->normal();
210 8839738 : const auto u_neighbor_normal = u_neighbor * _face_info->normal();
211 :
212 8839738 : const auto pr = computeAlphaAndOmega(u_elem_normal, u_neighbor_normal, c_elem, c_neighbor);
213 : const auto & alpha = pr.first;
214 : const auto & omega = pr.second;
215 :
216 8839738 : if (_eqn == "mass")
217 : return computeFaceFlux(
218 2337402 : alpha, omega, sup_vel_elem_normal, sup_vel_neighbor_normal, rho_elem, rho_neighbor);
219 6502336 : else if (_eqn == "momentum")
220 : {
221 4157684 : const auto rhou_elem = u_elem(_index) * rho_elem;
222 4157684 : const auto rhou_neighbor = u_neighbor(_index) * rho_neighbor;
223 4157684 : return computeFaceFlux(alpha,
224 : omega,
225 : sup_vel_elem_normal,
226 : sup_vel_neighbor_normal,
227 : rhou_elem,
228 : rhou_neighbor) +
229 8315368 : _face_info->normal()(_index) * (alpha * _eps_elem[_qp] * pressure_elem +
230 8315368 : (1 - alpha) * _eps_neighbor[_qp] * pressure_neighbor);
231 : }
232 2344652 : else if (_eqn == "energy")
233 : {
234 4674804 : const auto ht_elem = e_elem + 0.5 * u_elem * u_elem + pressure_elem / rho_elem;
235 : const auto ht_neighbor =
236 4674804 : e_neighbor + 0.5 * u_neighbor * u_neighbor + pressure_neighbor / rho_neighbor;
237 : const auto rho_ht_elem = rho_elem * ht_elem;
238 : const auto rho_ht_neighbor = rho_neighbor * ht_neighbor;
239 : return computeFaceFlux(
240 2337402 : alpha, omega, sup_vel_elem_normal, sup_vel_neighbor_normal, rho_ht_elem, rho_ht_neighbor);
241 : }
242 7250 : else if (_eqn == "scalar")
243 : {
244 : const auto scalar_elem = interpolate(*_limiter,
245 7250 : _scalar_elem[_qp],
246 7250 : _scalar_neighbor[_qp],
247 7250 : &(*_grad_scalar_elem)[_qp],
248 7250 : *_face_info,
249 7250 : true);
250 : const auto scalar_neighbor = interpolate(*_limiter,
251 7250 : _scalar_neighbor[_qp],
252 7250 : _scalar_elem[_qp],
253 7250 : &(*_grad_scalar_neighbor)[_qp],
254 7250 : *_face_info,
255 7250 : false);
256 : const auto rhos_elem = rho_elem * scalar_elem;
257 : const auto rhos_neighbor = rho_neighbor * scalar_neighbor;
258 : return computeFaceFlux(
259 7250 : alpha, omega, sup_vel_elem_normal, sup_vel_neighbor_normal, rhos_elem, rhos_neighbor);
260 : }
261 : else
262 0 : mooseError("Unrecognized enum type ", _eqn);
263 : }
|
