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 "FVFluxKernel.h"
11 :
12 : #include "MooseVariableFV.h"
13 : #include "SystemBase.h"
14 : #include "MooseMesh.h"
15 : #include "ADUtils.h"
16 : #include "RelationshipManager.h"
17 :
18 : #include "libmesh/elem.h"
19 : #include "libmesh/system.h"
20 :
21 : InputParameters
22 167395 : FVFluxKernel::validParams()
23 : {
24 167395 : InputParameters params = FVKernel::validParams();
25 167395 : params += TwoMaterialPropertyInterface::validParams();
26 167395 : params.registerSystemAttributeName("FVFluxKernel");
27 502185 : params.addParam<bool>("force_boundary_execution",
28 334790 : false,
29 : "Whether to force execution of this object on all external boundaries.");
30 502185 : params.addParam<std::vector<BoundaryName>>(
31 : "boundaries_to_force",
32 334790 : std::vector<BoundaryName>(),
33 : "The set of sidesets to force execution of this FVFluxKernel on. "
34 : "Setting force_boundary_execution to true is equivalent to listing all external "
35 : "mesh boundaries in this parameter.");
36 502185 : params.addParam<std::vector<BoundaryName>>(
37 : "boundaries_to_avoid",
38 334790 : std::vector<BoundaryName>(),
39 : "The set of sidesets to not execute this FVFluxKernel on. "
40 : "This takes precedence over force_boundary_execution to restrict to less external boundaries."
41 : " By default flux kernels are executed on all internal boundaries and Dirichlet boundary "
42 : "conditions.");
43 :
44 167395 : params.addParamNamesToGroup("force_boundary_execution boundaries_to_force boundaries_to_avoid",
45 : "Boundary execution modification");
46 167395 : return params;
47 0 : }
48 :
49 5205 : FVFluxKernel::FVFluxKernel(const InputParameters & params)
50 : : FVKernel(params),
51 : TwoMaterialPropertyInterface(this, blockIDs(), {}),
52 : NeighborMooseVariableInterface(
53 : this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
54 : NeighborCoupleableMooseVariableDependencyIntermediateInterface(
55 : this, false, false, /*is_fv=*/true),
56 10406 : _var(*mooseVariableFV()),
57 10402 : _u_elem(_var.adSln()),
58 5201 : _u_neighbor(_var.adSlnNeighbor()),
59 10406 : _force_boundary_execution(getParam<bool>("force_boundary_execution"))
60 : {
61 5201 : addMooseVariableDependency(&_var);
62 :
63 5201 : const auto & vec = getParam<std::vector<BoundaryName>>("boundaries_to_force");
64 5635 : for (const auto & name : vec)
65 434 : _boundaries_to_force.insert(_mesh.getBoundaryID(name));
66 :
67 5201 : const auto & avoid_vec = getParam<std::vector<BoundaryName>>("boundaries_to_avoid");
68 5227 : for (const auto & name : avoid_vec)
69 : {
70 30 : const auto bid = _mesh.getBoundaryID(name);
71 30 : _boundaries_to_avoid.insert(bid);
72 30 : if (_boundaries_to_force.find(bid) != _boundaries_to_force.end())
73 4 : paramError(
74 : "boundaries_to_avoid",
75 : "A boundary may not be specified in both boundaries_to_avoid and boundaries_to_force");
76 : }
77 5197 : }
78 :
79 : bool
80 29991038 : FVFluxKernel::onBoundary(const FaceInfo & fi) const
81 : {
82 29991038 : return Moose::FV::onBoundary(*this, fi);
83 : }
84 :
85 : // Note the lack of quadrature point loops in the residual/jacobian compute
86 : // functions. This is because finite volumes currently only works with
87 : // constant monomial elements. We only have one quadrature point regardless of
88 : // problem dimension and just multiply by the face area.
89 :
90 : bool
91 29991478 : FVFluxKernel::skipForBoundary(const FaceInfo & fi) const
92 : {
93 : // Boundaries to avoid come first, since they are always obeyed
94 29991478 : if (avoidBoundary(fi))
95 440 : return true;
96 :
97 : // We get this to check if we are on a kernel boundary or not
98 29991038 : const bool on_boundary = onBoundary(fi);
99 :
100 : // We are either on a kernel boundary or on an internal sideset
101 : // which is handled as a boundary
102 29991038 : if (on_boundary || !fi.boundaryIDs().empty())
103 : {
104 : // Blanket forcing on boundary
105 2170755 : if (_force_boundary_execution)
106 3400 : return false;
107 :
108 : // Selected boundaries to force
109 2170222 : for (const auto bnd_to_force : _boundaries_to_force)
110 4854 : if (fi.boundaryIDs().count(bnd_to_force))
111 1987 : return false;
112 :
113 : // If we have a flux boundary on this face, we skip. This
114 : // should be relatively easy to check with the cached maps.
115 2165368 : if (_var.getFluxBCs(fi).first)
116 710048 : return true;
117 :
118 : // If we have a dirichlet BC, we are not skipping
119 1455320 : if (_var.getDirichletBC(fi).first)
120 1084024 : return false;
121 : }
122 :
123 : // The last question is: are we on the inside or on the outside? If we are on an internal
124 : // face we dont skip, otherwise we assume a natural BC and skip
125 28191579 : return on_boundary;
126 : }
127 :
128 : void
129 22721641 : FVFluxKernel::computeResidual(const FaceInfo & fi)
130 : {
131 22721641 : if (skipForBoundary(fi))
132 794854 : return;
133 :
134 21926787 : _face_info = &fi;
135 21926787 : _normal = fi.normal();
136 21926787 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
137 21926787 : auto r = fi.faceArea() * fi.faceCoord() * MetaPhysicL::raw_value(computeQpResidual());
138 :
139 : // residual contributions for a flux kernel go to both neighboring faces.
140 : // They are equal in magnitude but opposite in direction due to the outward
141 : // facing unit normals of the face for each neighboring elements being
142 : // oriented oppositely. We calculate the residual contribution once using
143 : // the lower-id-elem-oriented _normal and just use the resulting residual's
144 : // negative for the contribution to the neighbor element.
145 :
146 : // The fancy face type if condition checks here are because we might
147 : // currently be running on a face for which this kernel's variable is only
148 : // defined on one side. If this is the case, we need to only calculate+add
149 : // the residual contribution if there is a dirichlet bc for the active
150 : // face+variable. We always need to add the residual contribution when the
151 : // variable is defined on both sides of the face. If the variable is only
152 : // defined on one side and there is NOT a dirichlet BC, then there is either
153 : // a flux BC or a natural BC - in either of those cases we don't want to add
154 : // any residual contributions from regular flux kernels.
155 21926787 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
156 21150593 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
157 : {
158 : // residual contribution of this kernel to the elem element
159 21921838 : prepareVectorTag(_assembly, _var.number());
160 21921838 : _local_re(0) = r;
161 21921838 : accumulateTaggedLocalResidual();
162 : }
163 21926787 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
164 21921838 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
165 : {
166 : // residual contribution of this kernel to the neighbor element
167 21150593 : prepareVectorTagNeighbor(_assembly, _var.number());
168 21150593 : _local_re(0) = -r;
169 21150593 : accumulateTaggedLocalResidual();
170 : }
171 : }
172 :
173 : void
174 7269837 : FVFluxKernel::computeJacobian(const FaceInfo & fi)
175 : {
176 7269837 : if (skipForBoundary(fi))
177 282405 : return;
178 :
179 6987432 : _face_info = &fi;
180 6987432 : _normal = fi.normal();
181 6987432 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
182 6987432 : const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
183 :
184 : // The fancy face type if condition checks here are because we might
185 : // currently be running on a face for which this kernel's variable is only
186 : // defined on one side. If this is the case, we need to only calculate+add
187 : // the jacobian contribution if there is a dirichlet bc for the active
188 : // face+variable. We always need to add the jacobian contribution when the
189 : // variable is defined on both sides of the face. If the variable is only
190 : // defined on one side and there is NOT a dirichlet BC, then there is either
191 : // a flux BC or a natural BC - in either of those cases we don't want to add
192 : // any jacobian contributions from regular flux kernels.
193 6987432 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
194 6679845 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
195 : {
196 : mooseAssert(_var.dofIndices().size() == 1, "We're currently built to use CONSTANT MONOMIALS");
197 :
198 27948876 : addResidualsAndJacobian(
199 13974438 : _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
200 : }
201 :
202 6987432 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
203 6987219 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
204 : {
205 : mooseAssert((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) ==
206 : (_var.dofIndices().size() == 0),
207 : "If the variable is only defined on the neighbor hand side of the face, then that "
208 : "means it should have no dof indices on the elem element. Conversely if "
209 : "the variable is defined on both sides of the face, then it should have a non-zero "
210 : "number of degrees of freedom on the elem element");
211 :
212 : // We switch the sign for the neighbor residual
213 6679845 : ADReal neighbor_r = -r;
214 :
215 : mooseAssert(_var.dofIndicesNeighbor().size() == 1,
216 : "We're currently built to use CONSTANT MONOMIALS");
217 :
218 26719380 : addResidualsAndJacobian(_assembly,
219 6679845 : std::array<ADReal, 1>{{neighbor_r}},
220 6679845 : _var.dofIndicesNeighbor(),
221 6679845 : _var.scalingFactor());
222 6679845 : }
223 20654496 : }
224 :
225 : void
226 156596 : FVFluxKernel::computeResidualAndJacobian(const FaceInfo & fi)
227 : {
228 156596 : computeJacobian(fi);
229 156596 : }
230 :
231 : ADReal
232 10686490 : FVFluxKernel::gradUDotNormal(const Moose::StateArg & time, const bool correct_skewness) const
233 : {
234 : mooseAssert(_face_info, "the face info should be non-null");
235 10686490 : return Moose::FV::gradUDotNormal(*_face_info, _var, time, correct_skewness);
236 : }
237 :
238 : Moose::ElemArg
239 10306759 : FVFluxKernel::elemArg(const bool correct_skewness) const
240 : {
241 : mooseAssert(_face_info, "the face info should be non-null");
242 10306759 : return {_face_info->elemPtr(), correct_skewness};
243 : }
244 :
245 : Moose::ElemArg
246 10306759 : FVFluxKernel::neighborArg(const bool correct_skewness) const
247 : {
248 : mooseAssert(_face_info, "the face info should be non-null");
249 10306759 : return {_face_info->neighborPtr(), correct_skewness};
250 : }
251 :
252 : Moose::FaceArg
253 425657 : FVFluxKernel::singleSidedFaceArg(const FaceInfo * fi,
254 : const Moose::FV::LimiterType limiter_type,
255 : const bool correct_skewness,
256 : const Moose::StateArg * state_limiter) const
257 : {
258 425657 : if (!fi)
259 425657 : fi = _face_info;
260 :
261 425657 : return makeFace(*fi, limiter_type, true, correct_skewness, state_limiter);
262 : }
263 :
264 : bool
265 29991478 : FVFluxKernel::avoidBoundary(const FaceInfo & fi) const
266 : {
267 32309894 : for (const auto bnd_id : fi.boundaryIDs())
268 2318856 : if (_boundaries_to_avoid.count(bnd_id))
269 440 : return true;
270 29991038 : return false;
271 : }
272 :
273 : void
274 0 : FVFluxKernel::computeResidual()
275 : {
276 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
277 : }
278 :
279 : void
280 0 : FVFluxKernel::computeJacobian()
281 : {
282 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
283 : }
284 :
285 : void
286 0 : FVFluxKernel::computeResidualAndJacobian()
287 : {
288 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
289 : }
290 :
291 : bool
292 29887950 : FVFluxKernel::hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const
293 : {
294 29887950 : if (fi_elem_side)
295 14943975 : return hasBlocks(fi.elem().subdomain_id());
296 : else
297 14943975 : return fi.neighborPtr() && hasBlocks(fi.neighbor().subdomain_id());
298 : }
|
