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 31666 : INSFVMomentumDiffusion::validParams()
24 : {
25 31666 : auto params = INSFVFluxKernel::validParams();
26 31666 : params += FVDiffusionInterpolationInterface::validParams();
27 31666 : params.addRequiredParam<MooseFunctorName>(NS::mu, "The viscosity");
28 31666 : params.addClassDescription(
29 : "Implements the Laplace form of the viscous stress in the Navier-Stokes equation.");
30 :
31 63332 : MooseEnum coeff_interp_method("average harmonic", "harmonic");
32 63332 : params.addParam<MooseEnum>("mu_interp_method",
33 : coeff_interp_method,
34 : "Switch that can select face interpolation method for the viscosity.");
35 31666 : 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 63332 : params.addRelationshipManager(
40 : "ElementSideNeighborLayers",
41 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
42 : Moose::RelationshipManagerType::COUPLING,
43 : [](const InputParameters & obj_params, InputParameters & rm_params)
44 36782 : { FVRelationshipManagerInterface::setRMParamsDiffusion(obj_params, rm_params, 3); });
45 :
46 63332 : params.addParam<bool>(
47 : "complete_expansion",
48 63332 : false,
49 : "Boolean parameter to use complete momentum expansion is the diffusion term.");
50 63332 : params.addParam<MooseFunctorName>("u", "The velocity in the x direction.");
51 63332 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
52 63332 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
53 63332 : params.addParam<bool>(
54 63332 : "limit_interpolation", false, "Flag to limit interpolation to positive values.");
55 63332 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
56 63332 : params.addParamNamesToGroup("newton_solve", "Advanced");
57 31666 : return params;
58 31666 : }
59 :
60 17593 : INSFVMomentumDiffusion::INSFVMomentumDiffusion(const InputParameters & params)
61 : : INSFVFluxKernel(params),
62 : SolutionInvalidInterface(this),
63 : FVDiffusionInterpolationInterface(params),
64 35186 : _mu(getFunctor<ADReal>(NS::mu)),
65 17593 : _mu_interp_method(
66 35186 : Moose::FV::selectInterpolationMethod(getParam<MooseEnum>("mu_interp_method"))),
67 21197 : _u_var(params.isParamValid("u") ? &getFunctor<ADReal>("u") : nullptr),
68 21197 : _v_var(params.isParamValid("v") ? &getFunctor<ADReal>("v") : nullptr),
69 17593 : _w_var(params.isParamValid("w") ? &getFunctor<ADReal>("w") : nullptr),
70 35186 : _complete_expansion(getParam<bool>("complete_expansion")),
71 35186 : _limit_interpolation(getParam<bool>("limit_interpolation")),
72 17593 : _dim(_subproblem.mesh().dimension()),
73 70372 : _newton_solve(getParam<bool>("newton_solve"))
74 : {
75 17593 : if (_complete_expansion && !_u_var)
76 0 : paramError("u", "The u velocity must be defined when 'complete_expansion=true'.");
77 :
78 17593 : if (_complete_expansion && _dim >= 2 && !_v_var)
79 0 : paramError("v",
80 : "The v velocity must be defined when 'complete_expansion=true'"
81 : "and problem dimension is larger or equal to 2.");
82 :
83 17593 : if (_complete_expansion && _dim >= 3 && !_w_var)
84 0 : paramError("w",
85 : "The w velocity must be defined when 'complete_expansion=true'"
86 : "and problem dimension is larger or equal to three.");
87 17593 : }
88 :
89 : ADReal
90 68129875 : INSFVMomentumDiffusion::computeStrongResidual(const bool populate_a_coeffs)
91 : {
92 68129875 : const Moose::StateArg state = determineState();
93 68129875 : const auto dudn = gradUDotNormal(state, _correct_skewness);
94 : ADReal face_mu;
95 :
96 68129875 : if (onBoundary(*_face_info))
97 6214730 : face_mu = _mu(makeCDFace(*_face_info), state);
98 : else
99 61915145 : Moose::FV::interpolate(_mu_interp_method,
100 : face_mu,
101 61915145 : _mu(elemArg(), state),
102 123830290 : _mu(neighborArg(), state),
103 61915145 : *_face_info,
104 : true);
105 :
106 : // Protecting from negative viscosity at interpolation
107 : // to preserve convergence
108 68129875 : if (face_mu < 0.0)
109 : {
110 216 : if (!(_limit_interpolation))
111 : {
112 108 : mooseDoOnce(mooseWarning(
113 : "Negative face viscosity has been encountered. Value ",
114 : raw_value(face_mu),
115 : " at ",
116 : _face_info->faceCentroid(),
117 : " limiting it to 0!\nFurther warnings for this issue will be silenced, but the "
118 : "occurrences will be recorded through the solution invalidity interface."));
119 117 : flagInvalidSolution("Negative face dynamic viscosity has been encountered.");
120 : }
121 : // Keep face_mu here for sparsity pattern detection
122 432 : face_mu = 0 * face_mu;
123 : }
124 :
125 68129875 : if (populate_a_coeffs)
126 : {
127 56621403 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
128 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
129 : {
130 56620593 : const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
131 : // A gradient is a linear combination of degrees of freedom so it's safe to straight-up index
132 : // into the derivatives vector at the dof we care about
133 56620593 : _ae = dudn.derivatives()[dof_number];
134 56620593 : _ae *= -face_mu;
135 : }
136 56621403 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
137 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
138 : {
139 52764394 : const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
140 52764394 : _an = dudn.derivatives()[dof_number];
141 52764394 : _an *= face_mu;
142 : }
143 : }
144 :
145 68129875 : ADReal dudn_transpose = 0.0;
146 68129875 : if (_complete_expansion)
147 : {
148 : // Computing the gradient from coupled variables
149 : // Normally, we can do this with `_var.gradient(face, state)` but we will need the transpose
150 : // gradient. So, we compute all at once
151 : Moose::FaceArg face;
152 2642440 : if (onBoundary(*_face_info))
153 399840 : face = singleSidedFaceArg();
154 : else
155 2242600 : face = makeCDFace(*_face_info, _correct_skewness);
156 :
157 : ADRealTensorValue gradient;
158 2642440 : if (_dim == 1)
159 : {
160 0 : const auto & grad_u = _u_var->gradient(face, state);
161 0 : gradient = ADRealTensorValue(grad_u, ADRealVectorValue(0, 0, 0), ADRealVectorValue(0, 0, 0));
162 : }
163 2642440 : else if (_dim == 2)
164 : {
165 2642440 : const auto & grad_u = _u_var->gradient(face, state);
166 2642440 : const auto & grad_v = _v_var->gradient(face, state);
167 5284880 : gradient = ADRealTensorValue(grad_u, grad_v, ADRealVectorValue(0, 0, 0));
168 : }
169 : else // if (_dim == 3)
170 : {
171 0 : const auto & grad_u = _u_var->gradient(face, state);
172 0 : const auto & grad_v = _v_var->gradient(face, state);
173 0 : const auto & grad_w = _w_var->gradient(face, state);
174 0 : gradient = ADRealTensorValue(grad_u, grad_v, grad_w);
175 : }
176 :
177 : // Getting transpose of the gradient matrix
178 2642440 : const auto gradient_transpose = gradient.transpose();
179 :
180 2642440 : dudn_transpose += gradient_transpose.row(_index) * _face_info->normal();
181 : }
182 :
183 136259750 : return -face_mu * (dudn + dudn_transpose);
184 : }
185 :
186 : void
187 96478312 : INSFVMomentumDiffusion::gatherRCData(const FaceInfo & fi)
188 : {
189 96478312 : if (skipForBoundary(fi))
190 : return;
191 :
192 94279455 : _face_info = &fi;
193 94279455 : _normal = fi.normal();
194 94279455 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
195 :
196 188558910 : addResidualAndJacobian(computeStrongResidual(true) * (fi.faceArea() * fi.faceCoord()));
197 :
198 94279455 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
199 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
200 188557050 : _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
201 94279455 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
202 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
203 267208002 : _rc_uo.addToA(fi.neighborPtr(), _index, _an * (fi.faceArea() * fi.faceCoord()));
204 : }
205 :
206 : ADReal
207 13803892 : INSFVMomentumDiffusion::computeSegregatedContribution()
208 : {
209 13803892 : return computeStrongResidual(false);
210 : }
|
