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 "INSFVTurbulentDiffusion.h"
11 : #include "NavierStokesMethods.h"
12 :
13 : registerMooseObject("NavierStokesApp", INSFVTurbulentDiffusion);
14 :
15 : InputParameters
16 1828 : INSFVTurbulentDiffusion::validParams()
17 : {
18 1828 : InputParameters params = FVDiffusion::validParams();
19 1828 : params.addClassDescription(
20 : "Computes residual for the turbulent scaled diffusion operator for finite volume method.");
21 3656 : params.addParam<MooseFunctorName>(
22 3656 : "scaling_coef", 1.0, "Scaling factor to divide the diffusion coefficient with");
23 3656 : params.addParam<std::vector<BoundaryName>>(
24 : "walls", {}, "Boundaries that correspond to solid walls.");
25 :
26 1828 : return params;
27 0 : }
28 :
29 1028 : INSFVTurbulentDiffusion::INSFVTurbulentDiffusion(const InputParameters & params)
30 : : FVDiffusion(params),
31 1028 : _scaling_coef(getFunctor<ADReal>("scaling_coef")),
32 2056 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
33 2056 : _preserve_sparsity_pattern(_fe_problem.preserveMatrixSparsityPattern())
34 : {
35 1028 : }
36 :
37 : void
38 3716 : INSFVTurbulentDiffusion::initialSetup()
39 : {
40 3716 : FVDiffusion::initialSetup();
41 7432 : NS::getWallBoundedElements(
42 3716 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
43 3716 : }
44 :
45 : ADReal
46 2700672 : INSFVTurbulentDiffusion::computeQpResidual()
47 : {
48 : using namespace Moose::FV;
49 2700672 : const auto state = determineState();
50 :
51 2700672 : const auto dudn = gradUDotNormal(state, _correct_skewness);
52 : ADReal coeff;
53 : ADReal scaling_coef;
54 :
55 : // If we are on internal faces, we interpolate the diffusivity as usual
56 2700672 : if (_var.isInternalFace(*_face_info))
57 : {
58 : // If the diffusion coefficients are zero, then we can early return 0 (and avoid warnings if we
59 : // have a harmonic interpolation)
60 2489936 : const auto coeff_elem = _coeff(elemArg(), state);
61 2489936 : const auto coeff_neighbor = _coeff(neighborArg(), state);
62 2489936 : if (!coeff_elem.value() && !coeff_neighbor.value())
63 : {
64 5280 : if (!_preserve_sparsity_pattern)
65 0 : return 0;
66 : else
67 : // we return 0 but preserve the sparsity pattern of the Jacobian for Newton's method
68 5280 : return 0 * (coeff_elem + coeff_neighbor) *
69 10560 : (_scaling_coef(elemArg(), state) + _scaling_coef(neighborArg(), state)) * dudn;
70 : }
71 :
72 2484656 : interpolate(_coeff_interp_method, coeff, coeff_elem, coeff_neighbor, *_face_info, true);
73 2484656 : interpolate(_coeff_interp_method,
74 : scaling_coef,
75 2484656 : _scaling_coef(elemArg(), state),
76 4969312 : _scaling_coef(neighborArg(), state),
77 2484656 : *_face_info,
78 : true);
79 : }
80 : // Else we just use the boundary values (which depend on how the diffusion
81 : // coefficient is constructed)
82 : else
83 : {
84 210736 : const auto face = singleSidedFaceArg();
85 210736 : coeff = _coeff(face, state);
86 210736 : scaling_coef = _scaling_coef(face, state);
87 : }
88 :
89 5390784 : return -1 * coeff / scaling_coef * dudn;
90 : }
91 :
92 : void
93 889184 : INSFVTurbulentDiffusion::computeResidual(const FaceInfo & fi)
94 : {
95 889184 : if (skipForBoundary(fi))
96 89216 : return;
97 :
98 799968 : _face_info = &fi;
99 799968 : _normal = fi.normal();
100 799968 : _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
101 1599936 : auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
102 :
103 799968 : const Elem * elem = fi.elemPtr();
104 799968 : const Elem * neighbor = fi.neighborPtr();
105 : const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
106 : const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
107 :
108 799968 : if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
109 799968 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
110 : (!bounded_elem))
111 : {
112 : // residual contribution of this kernel to the elem element
113 742144 : prepareVectorTag(_assembly, _var.number());
114 742144 : _local_re(0) = r;
115 742144 : accumulateTaggedLocalResidual();
116 : }
117 799968 : if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
118 740128 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
119 : (!bounded_neigh))
120 : {
121 : // residual contribution of this kernel to the neighbor element
122 682304 : prepareVectorTagNeighbor(_assembly, _var.number());
123 682304 : _local_re(0) = -r;
124 682304 : accumulateTaggedLocalResidual();
125 : }
126 : }
127 :
128 : void
129 2242448 : INSFVTurbulentDiffusion::computeJacobian(const FaceInfo & fi)
130 : {
131 2242448 : if (skipForBoundary(fi))
132 341744 : return;
133 :
134 1900704 : _face_info = &fi;
135 1900704 : _normal = fi.normal();
136 1900704 : _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
137 3801408 : const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
138 :
139 1900704 : const Elem * elem = fi.elemPtr();
140 1900704 : const Elem * neighbor = fi.neighborPtr();
141 : const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
142 : const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
143 :
144 1900704 : if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
145 1900704 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
146 : (!bounded_elem))
147 : {
148 : mooseAssert(_var.dofIndices().size() == 1, "We're currently built to use CONSTANT MONOMIALS");
149 :
150 1714568 : addResidualsAndJacobian(
151 3429136 : _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
152 : }
153 :
154 1900704 : if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
155 1749808 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
156 : (!bounded_neigh))
157 : {
158 : mooseAssert((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) ==
159 : (_var.dofIndices().size() == 0),
160 : "If the variable is only defined on the neighbor hand side of the face, then that "
161 : "means it should have no dof indices on the elem element. Conversely if "
162 : "the variable is defined on both sides of the face, then it should have a non-zero "
163 : "number of degrees of freedom on the elem element");
164 :
165 : // We switch the sign for the neighbor residual
166 1576288 : ADReal neighbor_r = -r;
167 :
168 : mooseAssert(_var.dofIndicesNeighbor().size() == 1,
169 : "We're currently built to use CONSTANT MONOMIALS");
170 :
171 1576288 : addResidualsAndJacobian(_assembly,
172 3152576 : std::array<ADReal, 1>{{neighbor_r}},
173 : _var.dofIndicesNeighbor(),
174 1576288 : _var.scalingFactor());
175 : }
176 : }
|
