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 "INSFVMomentumDiffusion.h"
11 : #include "INSFVRhieChowInterpolator.h"
12 : #include "NS.h"
13 : #include "NavierStokesMethods.h"
14 : #include "SystemBase.h"
15 : #include "NonlinearSystemBase.h"
16 : #include "RelationshipManager.h"
17 : #include "Factory.h"
18 : #include "libmesh/nonlinear_solver.h"
19 :
20 : registerMooseObject("NavierStokesApp", INSFVMomentumDiffusion);
21 :
22 : InputParameters
23 18369 : INSFVMomentumDiffusion::validParams()
24 : {
25 18369 : auto params = INSFVFluxKernel::validParams();
26 18369 : params += FVDiffusionInterpolationInterface::validParams();
27 18369 : params.addRequiredParam<MooseFunctorName>(NS::mu, "The viscosity");
28 18369 : params.addClassDescription(
29 : "Implements the Laplace form of the viscous stress in the Navier-Stokes equation.");
30 :
31 36738 : MooseEnum coeff_interp_method("average harmonic", "harmonic");
32 36738 : params.addParam<MooseEnum>("mu_interp_method",
33 : coeff_interp_method,
34 : "Switch that can select face interpolation method for the viscosity.");
35 18369 : params.set<unsigned short>("ghost_layers") = 2;
36 :
37 : // We add the relationship manager here, this will select the right number of
38 : // ghosting layers depending on the chosen interpolation method
39 36738 : params.addRelationshipManager(
40 : "ElementSideNeighborLayers",
41 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
42 : Moose::RelationshipManagerType::COUPLING,
43 : [](const InputParameters & obj_params, InputParameters & rm_params)
44 21845 : { FVRelationshipManagerInterface::setRMParamsDiffusion(obj_params, rm_params, 3); });
45 :
46 36738 : params.addParam<bool>(
47 : "complete_expansion",
48 36738 : false,
49 : "Boolean parameter to use complete momentum expansion is the diffusion term.");
50 36738 : params.addParam<bool>("include_isotropic_viscous_stress",
51 36738 : false,
52 : "Add the -(2/3) mu div(u) I term (requires specifying the velocity "
53 : "components via 'u', 'v', 'w'). Only meaningful for "
54 : "weakly-compressible formulations.");
55 36738 : params.addParam<MooseFunctorName>("u", "The velocity in the x direction.");
56 36738 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
57 36738 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
58 36738 : params.addParam<bool>(
59 36738 : "limit_interpolation", false, "Flag to limit interpolation to positive values.");
60 36738 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
61 36738 : params.addParamNamesToGroup("newton_solve", "Advanced");
62 18369 : return params;
63 18369 : }
64 :
65 10021 : INSFVMomentumDiffusion::INSFVMomentumDiffusion(const InputParameters & params)
66 : : INSFVFluxKernel(params),
67 : FVDiffusionInterpolationInterface(params),
68 20042 : _mu(getFunctor<ADReal>(NS::mu)),
69 10021 : _mu_interp_method(
70 20042 : Moose::FV::selectInterpolationMethod(getParam<MooseEnum>("mu_interp_method"))),
71 12057 : _u_var(params.isParamValid("u") ? &getFunctor<ADReal>("u") : nullptr),
72 12057 : _v_var(params.isParamValid("v") ? &getFunctor<ADReal>("v") : nullptr),
73 10021 : _w_var(params.isParamValid("w") ? &getFunctor<ADReal>("w") : nullptr),
74 20042 : _complete_expansion(getParam<bool>("complete_expansion")),
75 20042 : _include_isotropic_viscous_stress(getParam<bool>("include_isotropic_viscous_stress")),
76 20042 : _limit_interpolation(getParam<bool>("limit_interpolation")),
77 10021 : _dim(_subproblem.mesh().dimension()),
78 40084 : _newton_solve(getParam<bool>("newton_solve"))
79 : {
80 10021 : if (_complete_expansion && !_u_var)
81 0 : paramError("u", "The u velocity must be defined when 'complete_expansion=true'.");
82 :
83 10021 : if (_complete_expansion && _dim >= 2 && !_v_var)
84 0 : paramError("v",
85 : "The v velocity must be defined when 'complete_expansion=true'"
86 : "and problem dimension is larger or equal to 2.");
87 :
88 10021 : if (_complete_expansion && _dim >= 3 && !_w_var)
89 0 : paramError("w",
90 : "The w velocity must be defined when 'complete_expansion=true'"
91 : "and problem dimension is larger or equal to three.");
92 :
93 10021 : if (_include_isotropic_viscous_stress)
94 : {
95 0 : if (!_complete_expansion)
96 0 : paramError("include_isotropic_viscous_stress",
97 : "Complete expansion needs to be enabled to use the isotropic viscous stress!");
98 :
99 0 : if (!_u_var)
100 0 : paramError("include_isotropic_viscous_stress",
101 : "Velocity components must be provided to use the "
102 : "'include_isotropic_viscous_stress' option.");
103 0 : if (_dim >= 2 && !_v_var)
104 0 : paramError("include_isotropic_viscous_stress",
105 : "Velocity components must be provided to use the "
106 : "'include_isotropic_viscous_stress' option in dimensions >= 2.");
107 0 : if (_dim >= 3 && !_w_var)
108 0 : paramError("include_isotropic_viscous_stress",
109 : "Velocity components must be provided to use the "
110 : "'include_isotropic_viscous_stress' option in 3D.");
111 : }
112 10021 : }
113 :
114 : ADReal
115 48239701 : INSFVMomentumDiffusion::computeStrongResidual(const bool populate_a_coeffs)
116 : {
117 48239701 : const Moose::StateArg state = determineState();
118 48239701 : const auto dudn = gradUDotNormal(state, _correct_skewness);
119 : ADReal face_mu;
120 :
121 48239701 : if (onBoundary(*_face_info))
122 4435102 : face_mu = _mu(makeCDFace(*_face_info), state);
123 : else
124 43804599 : Moose::FV::interpolate(_mu_interp_method,
125 : face_mu,
126 43804599 : _mu(elemArg(), state),
127 87609198 : _mu(neighborArg(), state),
128 43804599 : *_face_info,
129 : true);
130 :
131 : // Protecting from negative viscosity at interpolation
132 : // to preserve convergence
133 48239701 : if (face_mu < 0.0)
134 : {
135 144 : if (!(_limit_interpolation))
136 : {
137 72 : mooseDoOnce(mooseWarning(
138 : "Negative face viscosity has been encountered. Value ",
139 : raw_value(face_mu),
140 : " at ",
141 : _face_info->faceCentroid(),
142 : " limiting it to 0!\nFurther warnings for this issue will be silenced, but the "
143 : "occurrences will be recorded through the solution invalidity interface."));
144 78 : flagInvalidSolution("Negative face dynamic viscosity has been encountered.");
145 : }
146 : // Keep face_mu here for sparsity pattern detection
147 288 : face_mu = 0 * face_mu;
148 : }
149 :
150 48239701 : if (populate_a_coeffs)
151 : {
152 40024607 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
153 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
154 : {
155 40024015 : const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
156 : // A gradient is a linear combination of degrees of freedom so it's safe to straight-up index
157 : // into the derivatives vector at the dof we care about
158 40024015 : _ae = dudn.derivatives()[dof_number];
159 40024015 : _ae *= -face_mu;
160 : }
161 40024607 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
162 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
163 : {
164 37202853 : const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
165 37202853 : _an = dudn.derivatives()[dof_number];
166 37202853 : _an *= face_mu;
167 : }
168 : }
169 :
170 48239701 : ADReal dudn_transpose = 0.0;
171 48239701 : ADReal divergence_term = 0.0;
172 48239701 : if (_complete_expansion)
173 : {
174 : // Computing the gradient from coupled variables
175 : // Normally, we can do this with `_var.gradient(face, state)` but we will need the transpose
176 : // gradient. So, we compute all at once
177 : Moose::FaceArg face;
178 1900408 : if (onBoundary(*_face_info))
179 283296 : face = singleSidedFaceArg();
180 : else
181 1617112 : face = makeCDFace(*_face_info, _correct_skewness);
182 :
183 : ADRealTensorValue gradient;
184 1900408 : if (_dim == 1)
185 : {
186 0 : const auto & grad_u = _u_var->gradient(face, state);
187 0 : gradient = ADRealTensorValue(grad_u, ADRealVectorValue(0, 0, 0), ADRealVectorValue(0, 0, 0));
188 : }
189 1900408 : else if (_dim == 2)
190 : {
191 1900408 : const auto & grad_u = _u_var->gradient(face, state);
192 1900408 : const auto & grad_v = _v_var->gradient(face, state);
193 3800816 : gradient = ADRealTensorValue(grad_u, grad_v, ADRealVectorValue(0, 0, 0));
194 : }
195 : else // if (_dim == 3)
196 : {
197 0 : const auto & grad_u = _u_var->gradient(face, state);
198 0 : const auto & grad_v = _v_var->gradient(face, state);
199 0 : const auto & grad_w = _w_var->gradient(face, state);
200 0 : gradient = ADRealTensorValue(grad_u, grad_v, grad_w);
201 : }
202 :
203 : // Getting transpose of the gradient matrix
204 1900408 : const auto gradient_transpose = gradient.transpose();
205 :
206 1900408 : dudn_transpose += gradient_transpose.row(_index) * _face_info->normal();
207 :
208 1900408 : if (_include_isotropic_viscous_stress)
209 : {
210 : libMesh::VectorValue<ADReal> velocity;
211 0 : velocity(0) = _u_var ? (*_u_var)(face, state) : ADReal(0);
212 0 : velocity(1) = _v_var ? (*_v_var)(face, state) : ADReal(0);
213 0 : velocity(2) = _w_var ? (*_w_var)(face, state) : ADReal(0);
214 :
215 0 : const auto coord_sys = _subproblem.getCoordSystem(_face_info->elem().subdomain_id());
216 : unsigned int rz_radial_coord =
217 0 : coord_sys == Moose::COORD_RZ ? _subproblem.getAxisymmetricRadialCoord() : 0;
218 : divergence_term =
219 0 : NS::divergence(gradient, velocity, face.getPoint(), coord_sys, rz_radial_coord);
220 : }
221 : }
222 :
223 48239701 : ADReal residual = -face_mu * (dudn + dudn_transpose);
224 48239701 : if (_complete_expansion && _include_isotropic_viscous_stress)
225 0 : residual += (2.0 / 3.0) * face_mu * divergence_term * _face_info->normal()(_index);
226 :
227 48239701 : return residual;
228 : }
229 :
230 : void
231 70355576 : INSFVMomentumDiffusion::gatherRCData(const FaceInfo & fi)
232 : {
233 70355576 : if (skipForBoundary(fi))
234 : return;
235 :
236 68708783 : _face_info = &fi;
237 68708783 : _normal = fi.normal();
238 68708783 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
239 :
240 137417566 : addResidualAndJacobian(computeStrongResidual(true) * (fi.faceArea() * fi.faceCoord()));
241 :
242 68708783 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
243 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
244 137416182 : _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
245 68708783 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
246 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
247 194820948 : _rc_uo.addToA(fi.neighborPtr(), _index, _an * (fi.faceArea() * fi.faceCoord()));
248 : }
249 :
250 : ADReal
251 10128594 : INSFVMomentumDiffusion::computeSegregatedContribution()
252 : {
253 10128594 : return computeStrongResidual(false);
254 : }
|
