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 3563 : PCNSFVKT::validParams()
21 : {
22 3563 : InputParameters params = FVFluxKernel::validParams();
23 3563 : params.addClassDescription("Computes the residual of advective term using finite volume method.");
24 3563 : params.addRequiredParam<UserObjectName>(NS::fluid, "Fluid userobject");
25 7126 : MooseEnum eqn("mass momentum energy scalar");
26 7126 : params.addRequiredParam<MooseEnum>("eqn", eqn, "The equation you're solving.");
27 7126 : MooseEnum momentum_component("x=0 y=1 z=2");
28 7126 : params.addParam<MooseEnum>("momentum_component",
29 : momentum_component,
30 : "The component of the momentum equation that this kernel applies to.");
31 7126 : params.addParam<MooseEnum>(
32 : "limiter", moose_limiter_type, "The limiter to apply during interpolation.");
33 3563 : params.set<unsigned short>("ghost_layers") = 2;
34 7126 : params.addParam<bool>(
35 : "knp_for_omega",
36 7126 : 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 3563 : return params;
40 3563 : }
41 :
42 1816 : PCNSFVKT::PCNSFVKT(const InputParameters & params)
43 : : FVFluxKernel(params),
44 3632 : _fluid(getUserObject<SinglePhaseFluidProperties>(NS::fluid)),
45 :
46 1816 : _sup_vel_x_elem(getADMaterialProperty<Real>(NS::superficial_velocity_x)),
47 1816 : _sup_vel_x_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_x)),
48 1816 : _grad_sup_vel_x_elem(
49 1816 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_x))),
50 1816 : _grad_sup_vel_x_neighbor(
51 1816 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_x))),
52 1816 : _sup_vel_y_elem(getADMaterialProperty<Real>(NS::superficial_velocity_y)),
53 1816 : _sup_vel_y_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_y)),
54 1816 : _grad_sup_vel_y_elem(
55 1816 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_y))),
56 1816 : _grad_sup_vel_y_neighbor(
57 1816 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_y))),
58 1816 : _sup_vel_z_elem(getADMaterialProperty<Real>(NS::superficial_velocity_z)),
59 1816 : _sup_vel_z_neighbor(getNeighborADMaterialProperty<Real>(NS::superficial_velocity_z)),
60 1816 : _grad_sup_vel_z_elem(
61 1816 : getADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_z))),
62 1816 : _grad_sup_vel_z_neighbor(
63 1816 : getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::superficial_velocity_z))),
64 :
65 1816 : _T_fluid_elem(getADMaterialProperty<Real>(NS::T_fluid)),
66 1816 : _T_fluid_neighbor(getNeighborADMaterialProperty<Real>(NS::T_fluid)),
67 1816 : _grad_T_fluid_elem(getADMaterialProperty<RealVectorValue>(NS::grad(NS::T_fluid))),
68 1816 : _grad_T_fluid_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::T_fluid))),
69 1816 : _pressure_elem(getADMaterialProperty<Real>(NS::pressure)),
70 1816 : _pressure_neighbor(getNeighborADMaterialProperty<Real>(NS::pressure)),
71 1816 : _grad_pressure_elem(getADMaterialProperty<RealVectorValue>(NS::grad(NS::pressure))),
72 1816 : _grad_pressure_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::grad(NS::pressure))),
73 1816 : _eps_elem(getMaterialProperty<Real>(NS::porosity)),
74 1816 : _eps_neighbor(getNeighborMaterialProperty<Real>(NS::porosity)),
75 3632 : _eqn(getParam<MooseEnum>("eqn")),
76 3632 : _index(getParam<MooseEnum>("momentum_component")),
77 1816 : _scalar_elem(_u_elem),
78 1816 : _scalar_neighbor(_u_neighbor),
79 1816 : _grad_scalar_elem((_eqn == "scalar") ? &_var.adGradSln() : nullptr),
80 1816 : _grad_scalar_neighbor((_eqn == "scalar") ? &_var.adGradSlnNeighbor() : nullptr),
81 3632 : _limiter(Limiter<ADReal>::build(LimiterType(int(getParam<MooseEnum>("limiter"))))),
82 5448 : _knp_for_omega(getParam<bool>("knp_for_omega"))
83 : {
84 3745 : 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 1816 : }
89 :
90 : std::pair<ADReal, ADReal>
91 11567329 : 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 11567329 : std::max({c_elem + u_elem_normal, c_neighbor + u_neighbor_normal, ADReal(0)});
100 : const auto psi_neighbor =
101 11567329 : std::max({c_elem - u_elem_normal, c_neighbor - u_neighbor_normal, ADReal(0)});
102 22944218 : auto alpha = _knp_for_omega ? psi_elem / (psi_elem + psi_neighbor) : ADReal(0.5);
103 57265325 : auto omega = _knp_for_omega ? alpha * (1 - alpha) * (psi_elem + psi_neighbor)
104 11948209 : : alpha * std::max(psi_elem, psi_neighbor);
105 :
106 : // Do this to avoid new nonzero mallocs
107 23134658 : const auto dummy_quant = 0 * (c_elem + u_elem_normal + c_neighbor + u_neighbor_normal);
108 :
109 11567329 : alpha += dummy_quant;
110 11567329 : omega += dummy_quant;
111 11567329 : return std::make_pair(std::move(alpha), std::move(omega));
112 23134658 : }
113 :
114 : ADReal
115 11567329 : 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 11567329 : return alpha * (sup_vel_elem_normal * adv_quant_elem) +
123 11567329 : (1 - alpha) * sup_vel_neighbor_normal * adv_quant_neighbor -
124 11567329 : omega * (adv_quant_neighbor - adv_quant_elem);
125 : }
126 :
127 : ADReal
128 11567329 : PCNSFVKT::computeQpResidual()
129 : {
130 : // Perform primitive interpolations
131 : const auto pressure_elem = interpolate(*_limiter,
132 11567329 : _pressure_elem[_qp],
133 11567329 : _pressure_neighbor[_qp],
134 11567329 : &_grad_pressure_elem[_qp],
135 11567329 : *_face_info,
136 11567329 : /*elem_is_up=*/true);
137 : const auto pressure_neighbor = interpolate(*_limiter,
138 11567329 : _pressure_neighbor[_qp],
139 11567329 : _pressure_elem[_qp],
140 11567329 : &_grad_pressure_neighbor[_qp],
141 11567329 : *_face_info,
142 11567329 : /*elem_is_up=*/false);
143 : const auto T_fluid_elem = interpolate(*_limiter,
144 11567329 : _T_fluid_elem[_qp],
145 11567329 : _T_fluid_neighbor[_qp],
146 11567329 : &_grad_T_fluid_elem[_qp],
147 11567329 : *_face_info,
148 11567329 : /*elem_is_up=*/true);
149 : const auto T_fluid_neighbor = interpolate(*_limiter,
150 11567329 : _T_fluid_neighbor[_qp],
151 11567329 : _T_fluid_elem[_qp],
152 11567329 : &_grad_T_fluid_neighbor[_qp],
153 11567329 : *_face_info,
154 11567329 : /*elem_is_up=*/false);
155 : const auto sup_vel_x_elem = interpolate(*_limiter,
156 11567329 : _sup_vel_x_elem[_qp],
157 11567329 : _sup_vel_x_neighbor[_qp],
158 11567329 : &_grad_sup_vel_x_elem[_qp],
159 11567329 : *_face_info,
160 11567329 : /*elem_is_up=*/true);
161 : const auto sup_vel_x_neighbor = interpolate(*_limiter,
162 11567329 : _sup_vel_x_neighbor[_qp],
163 11567329 : _sup_vel_x_elem[_qp],
164 11567329 : &_grad_sup_vel_x_neighbor[_qp],
165 11567329 : *_face_info,
166 11567329 : /*elem_is_up=*/false);
167 : const auto sup_vel_y_elem = interpolate(*_limiter,
168 11567329 : _sup_vel_y_elem[_qp],
169 11567329 : _sup_vel_y_neighbor[_qp],
170 11567329 : &_grad_sup_vel_y_elem[_qp],
171 11567329 : *_face_info,
172 11567329 : /*elem_is_up=*/true);
173 : const auto sup_vel_y_neighbor = interpolate(*_limiter,
174 11567329 : _sup_vel_y_neighbor[_qp],
175 11567329 : _sup_vel_y_elem[_qp],
176 11567329 : &_grad_sup_vel_y_neighbor[_qp],
177 11567329 : *_face_info,
178 11567329 : /*elem_is_up=*/false);
179 : const auto sup_vel_z_elem = interpolate(*_limiter,
180 11567329 : _sup_vel_z_elem[_qp],
181 11567329 : _sup_vel_z_neighbor[_qp],
182 11567329 : &_grad_sup_vel_z_elem[_qp],
183 11567329 : *_face_info,
184 11567329 : /*elem_is_up=*/true);
185 : const auto sup_vel_z_neighbor = interpolate(*_limiter,
186 11567329 : _sup_vel_z_neighbor[_qp],
187 11567329 : _sup_vel_z_elem[_qp],
188 11567329 : &_grad_sup_vel_z_neighbor[_qp],
189 11567329 : *_face_info,
190 11567329 : /*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 11567329 : const auto u_elem = sup_vel_elem / _eps_elem[_qp];
194 11567329 : const auto rho_elem = _fluid.rho_from_p_T(pressure_elem, T_fluid_elem);
195 11567329 : const auto specific_volume_elem = 1. / rho_elem;
196 11567329 : 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 11567329 : const auto u_neighbor = sup_vel_neighbor / _eps_neighbor[_qp];
200 11567329 : const auto rho_neighbor = _fluid.rho_from_p_T(pressure_neighbor, T_fluid_neighbor);
201 11567329 : const auto specific_volume_neighbor = 1. / rho_neighbor;
202 11567329 : const auto e_neighbor = _fluid.e_from_T_v(T_fluid_neighbor, specific_volume_neighbor);
203 :
204 11567329 : const auto c_elem = _fluid.c_from_v_e(specific_volume_elem, e_elem);
205 11567329 : const auto c_neighbor = _fluid.c_from_v_e(specific_volume_neighbor, e_neighbor);
206 :
207 11567329 : const auto sup_vel_elem_normal = sup_vel_elem * _face_info->normal();
208 11567329 : const auto sup_vel_neighbor_normal = sup_vel_neighbor * _face_info->normal();
209 11567329 : const auto u_elem_normal = u_elem * _face_info->normal();
210 11567329 : const auto u_neighbor_normal = u_neighbor * _face_info->normal();
211 :
212 11567329 : 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 11567329 : if (_eqn == "mass")
217 : return computeFaceFlux(
218 3085759 : alpha, omega, sup_vel_elem_normal, sup_vel_neighbor_normal, rho_elem, rho_neighbor);
219 8481570 : else if (_eqn == "momentum")
220 : {
221 5384211 : const auto rhou_elem = u_elem(_index) * rho_elem;
222 5384211 : const auto rhou_neighbor = u_neighbor(_index) * rho_neighbor;
223 5384211 : return computeFaceFlux(alpha,
224 : omega,
225 : sup_vel_elem_normal,
226 : sup_vel_neighbor_normal,
227 : rhou_elem,
228 : rhou_neighbor) +
229 10768422 : _face_info->normal()(_index) * (alpha * _eps_elem[_qp] * pressure_elem +
230 10768422 : (1 - alpha) * _eps_neighbor[_qp] * pressure_neighbor);
231 : }
232 3097359 : else if (_eqn == "energy")
233 : {
234 6171518 : const auto ht_elem = e_elem + 0.5 * u_elem * u_elem + pressure_elem / rho_elem;
235 : const auto ht_neighbor =
236 6171518 : 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 3085759 : alpha, omega, sup_vel_elem_normal, sup_vel_neighbor_normal, rho_ht_elem, rho_ht_neighbor);
241 : }
242 11600 : else if (_eqn == "scalar")
243 : {
244 : const auto scalar_elem = interpolate(*_limiter,
245 11600 : _scalar_elem[_qp],
246 11600 : _scalar_neighbor[_qp],
247 11600 : &(*_grad_scalar_elem)[_qp],
248 11600 : *_face_info,
249 11600 : true);
250 : const auto scalar_neighbor = interpolate(*_limiter,
251 11600 : _scalar_neighbor[_qp],
252 11600 : _scalar_elem[_qp],
253 11600 : &(*_grad_scalar_neighbor)[_qp],
254 11600 : *_face_info,
255 11600 : false);
256 : const auto rhos_elem = rho_elem * scalar_elem;
257 : const auto rhos_neighbor = rho_neighbor * scalar_neighbor;
258 : return computeFaceFlux(
259 11600 : 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 : }
|
