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 3776 : INSFVTurbulentDiffusion::validParams()
17 : {
18 3776 : InputParameters params = FVDiffusion::validParams();
19 3776 : params.addClassDescription(
20 : "Computes residual for the turbulent scaled diffusion operator for finite volume method.");
21 7552 : params.addParam<MooseFunctorName>(
22 7552 : "scaling_coef", 1.0, "Scaling factor to divide the diffusion coefficient with");
23 7552 : params.addParam<std::vector<BoundaryName>>(
24 : "walls", {}, "Boundaries that correspond to solid walls.");
25 :
26 3776 : return params;
27 0 : }
28 :
29 2164 : INSFVTurbulentDiffusion::INSFVTurbulentDiffusion(const InputParameters & params)
30 : : FVDiffusion(params),
31 2164 : _scaling_coef(getFunctor<ADReal>("scaling_coef")),
32 4328 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
33 4328 : _preserve_sparsity_pattern(_fe_problem.preserveMatrixSparsityPattern())
34 : {
35 2164 : }
36 :
37 : void
38 7988 : INSFVTurbulentDiffusion::initialSetup()
39 : {
40 7988 : FVDiffusion::initialSetup();
41 15976 : NS::getWallBoundedElements(
42 7988 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
43 7988 : }
44 :
45 : ADReal
46 3722016 : INSFVTurbulentDiffusion::computeQpResidual()
47 : {
48 : using namespace Moose::FV;
49 3722016 : const auto state = determineState();
50 :
51 3722016 : 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 3722016 : 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 3409632 : const auto coeff_elem = _coeff(elemArg(), state);
61 3409632 : const auto coeff_neighbor = _coeff(neighborArg(), state);
62 3409632 : if (!coeff_elem.value() && !coeff_neighbor.value())
63 : {
64 7920 : 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 7920 : return 0 * (coeff_elem + coeff_neighbor) *
69 15840 : (_scaling_coef(elemArg(), state) + _scaling_coef(neighborArg(), state)) * dudn;
70 : }
71 :
72 3401712 : interpolate(_coeff_interp_method, coeff, coeff_elem, coeff_neighbor, *_face_info, true);
73 3401712 : interpolate(_coeff_interp_method,
74 : scaling_coef,
75 3401712 : _scaling_coef(elemArg(), state),
76 6803424 : _scaling_coef(neighborArg(), state),
77 3401712 : *_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 312384 : const auto face = singleSidedFaceArg();
85 312384 : coeff = _coeff(face, state);
86 312384 : scaling_coef = _scaling_coef(face, state);
87 : }
88 :
89 7428192 : return -1 * coeff / scaling_coef * dudn;
90 : }
91 :
92 : void
93 997760 : INSFVTurbulentDiffusion::computeResidual(const FaceInfo & fi)
94 : {
95 997760 : if (skipForBoundary(fi))
96 82720 : return;
97 :
98 915040 : _face_info = &fi;
99 915040 : _normal = fi.normal();
100 915040 : _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
101 1830080 : auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
102 :
103 915040 : const Elem * elem = fi.elemPtr();
104 915040 : 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 915040 : if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
109 915040 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
110 : (!bounded_elem))
111 : {
112 : // residual contribution of this kernel to the elem element
113 861960 : prepareVectorTag(_assembly, _var.number());
114 861960 : _local_re(0) = r;
115 861960 : accumulateTaggedLocalResidual();
116 : }
117 915040 : if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
118 826880 : _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
119 : (!bounded_neigh))
120 : {
121 : // residual contribution of this kernel to the neighbor element
122 773800 : prepareVectorTagNeighbor(_assembly, _var.number());
123 773800 : _local_re(0) = -r;
124 773800 : accumulateTaggedLocalResidual();
125 : }
126 : }
127 :
128 : void
129 3313952 : INSFVTurbulentDiffusion::computeJacobian(const FaceInfo & fi)
130 : {
131 3313952 : if (skipForBoundary(fi))
132 506976 : return;
133 :
134 2806976 : _face_info = &fi;
135 2806976 : _normal = fi.normal();
136 2806976 : _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
137 5613952 : const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
138 :
139 2806976 : const Elem * elem = fi.elemPtr();
140 2806976 : 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 2806976 : if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
145 2806976 : _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 2531272 : addResidualsAndJacobian(
151 5062544 : _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
152 : }
153 :
154 2806976 : if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
155 2582752 : _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 2325832 : ADReal neighbor_r = -r;
167 :
168 : mooseAssert(_var.dofIndicesNeighbor().size() == 1,
169 : "We're currently built to use CONSTANT MONOMIALS");
170 :
171 2325832 : addResidualsAndJacobian(_assembly,
172 4651664 : std::array<ADReal, 1>{{neighbor_r}},
173 : _var.dofIndicesNeighbor(),
174 2325832 : _var.scalingFactor());
175 : }
176 : }
|
