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 31896 : INSFVMomentumDiffusion::validParams()
24 : {
25 31896 : auto params = INSFVFluxKernel::validParams();
26 31896 : params += FVDiffusionInterpolationInterface::validParams();
27 31896 : params.addRequiredParam<MooseFunctorName>(NS::mu, "The viscosity");
28 31896 : params.addClassDescription(
29 : "Implements the Laplace form of the viscous stress in the Navier-Stokes equation.");
30 :
31 63792 : MooseEnum coeff_interp_method("average harmonic", "harmonic");
32 63792 : params.addParam<MooseEnum>("mu_interp_method",
33 : coeff_interp_method,
34 : "Switch that can select face interpolation method for the viscosity.");
35 31896 : 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 63792 : 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 63792 : params.addParam<bool>(
47 : "complete_expansion",
48 63792 : false,
49 : "Boolean parameter to use complete momentum expansion is the diffusion term.");
50 63792 : params.addParam<MooseFunctorName>("u", "The velocity in the x direction.");
51 63792 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
52 63792 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
53 63792 : params.addParam<bool>(
54 63792 : "limit_interpolation", false, "Flag to limit interpolation to positive values.");
55 63792 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
56 63792 : params.addParamNamesToGroup("newton_solve", "Advanced");
57 31896 : return params;
58 31896 : }
59 :
60 17747 : INSFVMomentumDiffusion::INSFVMomentumDiffusion(const InputParameters & params)
61 : : INSFVFluxKernel(params),
62 : FVDiffusionInterpolationInterface(params),
63 35494 : _mu(getFunctor<ADReal>(NS::mu)),
64 17747 : _mu_interp_method(
65 35494 : Moose::FV::selectInterpolationMethod(getParam<MooseEnum>("mu_interp_method"))),
66 21351 : _u_var(params.isParamValid("u") ? &getFunctor<ADReal>("u") : nullptr),
67 21351 : _v_var(params.isParamValid("v") ? &getFunctor<ADReal>("v") : nullptr),
68 17747 : _w_var(params.isParamValid("w") ? &getFunctor<ADReal>("w") : nullptr),
69 35494 : _complete_expansion(getParam<bool>("complete_expansion")),
70 35494 : _limit_interpolation(getParam<bool>("limit_interpolation")),
71 17747 : _dim(_subproblem.mesh().dimension()),
72 70988 : _newton_solve(getParam<bool>("newton_solve"))
73 : {
74 17747 : if (_complete_expansion && !_u_var)
75 0 : paramError("u", "The u velocity must be defined when 'complete_expansion=true'.");
76 :
77 17747 : 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 17747 : 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 17747 : }
87 :
88 : ADReal
89 69497360 : INSFVMomentumDiffusion::computeStrongResidual(const bool populate_a_coeffs)
90 : {
91 69497360 : const Moose::StateArg state = determineState();
92 69497360 : const auto dudn = gradUDotNormal(state, _correct_skewness);
93 : ADReal face_mu;
94 :
95 69497360 : if (onBoundary(*_face_info))
96 6382708 : face_mu = _mu(makeCDFace(*_face_info), state);
97 : else
98 63114652 : Moose::FV::interpolate(_mu_interp_method,
99 : face_mu,
100 63114652 : _mu(elemArg(), state),
101 126229304 : _mu(neighborArg(), state),
102 63114652 : *_face_info,
103 : true);
104 :
105 : // Protecting from negative viscosity at interpolation
106 : // to preserve convergence
107 69497360 : 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 69497360 : if (populate_a_coeffs)
125 : {
126 57989096 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
127 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
128 : {
129 57988286 : 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 57988286 : _ae = dudn.derivatives()[dof_number];
133 57988286 : _ae *= -face_mu;
134 : }
135 57989096 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
136 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
137 : {
138 53964053 : const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
139 53964053 : _an = dudn.derivatives()[dof_number];
140 53964053 : _an *= face_mu;
141 : }
142 : }
143 :
144 69497360 : ADReal dudn_transpose = 0.0;
145 69497360 : 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 138994720 : return -face_mu * (dudn + dudn_transpose);
183 : }
184 :
185 : void
186 97866637 : INSFVMomentumDiffusion::gatherRCData(const FaceInfo & fi)
187 : {
188 97866637 : if (skipForBoundary(fi))
189 : return;
190 :
191 95647148 : _face_info = &fi;
192 95647148 : _normal = fi.normal();
193 95647148 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
194 :
195 191294296 : addResidualAndJacobian(computeStrongResidual(true) * (fi.faceArea() * fi.faceCoord()));
196 :
197 95647148 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
198 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
199 191292436 : _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
200 95647148 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
201 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
202 270806979 : _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 : }
|
