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