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 42302 : FVFluxKernel::validParams()
23 : {
24 42302 : InputParameters params = FVKernel::validParams();
25 42302 : params += TwoMaterialPropertyInterface::validParams();
26 84604 : params.registerSystemAttributeName("FVFluxKernel");
27 126906 : params.addParam<bool>("force_boundary_execution",
28 84604 : false,
29 : "Whether to force execution of this object on all external boundaries.");
30 126906 : params.addParam<std::vector<BoundaryName>>(
31 : "boundaries_to_force",
32 84604 : 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 126906 : params.addParam<std::vector<BoundaryName>>(
37 : "boundaries_to_avoid",
38 84604 : 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 126906 : params.addParamNamesToGroup("force_boundary_execution boundaries_to_force boundaries_to_avoid",
45 : "Boundary execution modification");
46 42302 : return params;
47 0 : }
48 :
49 4143 : 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 8283 : _var(*mooseVariableFV()),
57 12423 : _force_boundary_execution(getParam<bool>("force_boundary_execution"))
58 : {
59 4140 : addMooseVariableDependency(&_var);
60 :
61 8280 : const auto & vec = getParam<std::vector<BoundaryName>>("boundaries_to_force");
62 4409 : for (const auto & name : vec)
63 269 : _boundaries_to_force.insert(_mesh.getBoundaryID(name));
64 :
65 8280 : const auto & avoid_vec = getParam<std::vector<BoundaryName>>("boundaries_to_avoid");
66 4166 : for (const auto & name : avoid_vec)
67 : {
68 29 : const auto bid = _mesh.getBoundaryID(name);
69 29 : _boundaries_to_avoid.insert(bid);
70 29 : if (_boundaries_to_force.find(bid) != _boundaries_to_force.end())
71 6 : paramError(
72 : "boundaries_to_avoid",
73 : "A boundary may not be specified in both boundaries_to_avoid and boundaries_to_force");
74 : }
75 4137 : }
76 :
77 : bool
78 18633008 : FVFluxKernel::onBoundary(const FaceInfo & fi) const
79 : {
80 18633008 : return Moose::FV::onBoundary(*this, fi);
81 : }
82 :
83 : // Note the lack of quadrature point loops in the residual/jacobian compute
84 : // functions. This is because finite volumes currently only works with
85 : // constant monomial elements. We only have one quadrature point regardless of
86 : // problem dimension and just multiply by the face area.
87 :
88 : bool
89 18633448 : FVFluxKernel::skipForBoundary(const FaceInfo & fi) const
90 : {
91 : // Boundaries to avoid come first, since they are always obeyed
92 18633448 : if (avoidBoundary(fi))
93 440 : return true;
94 :
95 : // We get this to check if we are on a kernel boundary or not
96 18633008 : const bool on_boundary = onBoundary(fi);
97 :
98 : // We are either on a kernel boundary or on an internal sideset
99 : // which is handled as a boundary
100 18633008 : if (on_boundary || !fi.boundaryIDs().empty())
101 : {
102 : // Blanket forcing on boundary
103 1742005 : if (_force_boundary_execution)
104 3064 : return false;
105 :
106 : // Selected boundaries to force
107 1741068 : for (const auto bnd_to_force : _boundaries_to_force)
108 3374 : if (fi.boundaryIDs().count(bnd_to_force))
109 1247 : return false;
110 :
111 : // If we have a flux boundary on this face, we skip. This
112 : // should be relatively easy to check with the cached maps.
113 1737694 : if (_var.getFluxBCs(fi).first)
114 536078 : return true;
115 :
116 : // If we have a dirichlet BC, we are not skipping
117 1201616 : if (_var.getDirichletBC(fi).first)
118 836153 : return false;
119 : }
120 :
121 : // The last question is: are we on the inside or on the outside? If we are on an internal
122 : // face we dont skip, otherwise we assume a natural BC and skip
123 17256466 : return on_boundary;
124 : }
125 :
126 : void
127 13782050 : FVFluxKernel::computeResidual(const FaceInfo & fi)
128 : {
129 13782050 : if (skipForBoundary(fi))
130 641901 : return;
131 :
132 13140149 : _face_info = &fi;
133 13140149 : _normal = fi.normal();
134 13140149 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
135 13140149 : auto r = fi.faceArea() * fi.faceCoord() * MetaPhysicL::raw_value(computeQpResidual());
136 :
137 : // residual contributions for a flux kernel go to both neighboring faces.
138 : // They are equal in magnitude but opposite in direction due to the outward
139 : // facing unit normals of the face for each neighboring elements being
140 : // oriented oppositely. We calculate the residual contribution once using
141 : // the lower-id-elem-oriented _normal and just use the resulting residual's
142 : // negative for the contribution to the neighbor element.
143 :
144 : // The fancy face type if condition checks here are because we might
145 : // currently be running on a face for which this kernel's variable is only
146 : // defined on one side. If this is the case, we need to only calculate+add
147 : // the residual contribution if there is a dirichlet bc for the active
148 : // face+variable. We always need to add the residual contribution when the
149 : // variable is defined on both sides of the face. If the variable is only
150 : // defined on one side and there is NOT a dirichlet BC, then there is either
151 : // a flux BC or a natural BC - in either of those cases we don't want to add
152 : // any residual contributions from regular flux kernels.
153 13140149 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
154 12558131 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
155 : {
156 : // residual contribution of this kernel to the elem element
157 13135981 : prepareVectorTag(_assembly, _var.number());
158 13135981 : _local_re(0) = r;
159 13135981 : accumulateTaggedLocalResidual();
160 : }
161 13140149 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
162 13135981 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
163 : {
164 : // residual contribution of this kernel to the neighbor element
165 12558131 : prepareVectorTagNeighbor(_assembly, _var.number());
166 12558131 : _local_re(0) = -r;
167 12558131 : accumulateTaggedLocalResidual();
168 : }
169 : }
170 :
171 : void
172 4851398 : FVFluxKernel::computeJacobian(const FaceInfo & fi)
173 : {
174 4851398 : if (skipForBoundary(fi))
175 255489 : return;
176 :
177 4595909 : _face_info = &fi;
178 4595909 : _normal = fi.normal();
179 4595909 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
180 4595909 : const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
181 :
182 : // The fancy face type if condition checks here are because we might
183 : // currently be running on a face for which this kernel's variable is only
184 : // defined on one side. If this is the case, we need to only calculate+add
185 : // the jacobian contribution if there is a dirichlet bc for the active
186 : // face+variable. We always need to add the jacobian contribution when the
187 : // variable is defined on both sides of the face. If the variable is only
188 : // defined on one side and there is NOT a dirichlet BC, then there is either
189 : // a flux BC or a natural BC - in either of those cases we don't want to add
190 : // any jacobian contributions from regular flux kernels.
191 4595909 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
192 4342248 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
193 : {
194 : mooseAssert(_var.dofIndices().size() == 1, "We're currently built to use CONSTANT MONOMIALS");
195 :
196 18383040 : addResidualsAndJacobian(
197 4595760 : _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
198 : }
199 :
200 4595909 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
201 4595760 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
202 : {
203 : mooseAssert((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) ==
204 : (_var.dofIndices().size() == 0),
205 : "If the variable is only defined on the neighbor hand side of the face, then that "
206 : "means it should have no dof indices on the elem element. Conversely if "
207 : "the variable is defined on both sides of the face, then it should have a non-zero "
208 : "number of degrees of freedom on the elem element");
209 :
210 : // We switch the sign for the neighbor residual
211 4342248 : ADReal neighbor_r = -r;
212 :
213 : mooseAssert(_var.dofIndicesNeighbor().size() == 1,
214 : "We're currently built to use CONSTANT MONOMIALS");
215 :
216 17368992 : addResidualsAndJacobian(_assembly,
217 4342248 : std::array<ADReal, 1>{{neighbor_r}},
218 4342248 : _var.dofIndicesNeighbor(),
219 4342248 : _var.scalingFactor());
220 4342248 : }
221 4595909 : }
222 :
223 : void
224 157220 : FVFluxKernel::computeResidualAndJacobian(const FaceInfo & fi)
225 : {
226 157220 : computeJacobian(fi);
227 157220 : }
228 :
229 : ADReal
230 8896737 : FVFluxKernel::gradUDotNormal(const Moose::StateArg & time, const bool correct_skewness) const
231 : {
232 : mooseAssert(_face_info, "the face info should be non-null");
233 8896737 : return Moose::FV::gradUDotNormal(*_face_info, _var, time, correct_skewness);
234 : }
235 :
236 : Moose::ElemArg
237 8579469 : FVFluxKernel::elemArg(const bool correct_skewness) const
238 : {
239 : mooseAssert(_face_info, "the face info should be non-null");
240 8579469 : return {_face_info->elemPtr(), correct_skewness};
241 : }
242 :
243 : Moose::ElemArg
244 8579469 : FVFluxKernel::neighborArg(const bool correct_skewness) const
245 : {
246 : mooseAssert(_face_info, "the face info should be non-null");
247 8579469 : return {_face_info->neighborPtr(), correct_skewness};
248 : }
249 :
250 : Moose::FaceArg
251 356816 : FVFluxKernel::singleSidedFaceArg(const FaceInfo * fi,
252 : const Moose::FV::LimiterType limiter_type,
253 : const bool correct_skewness,
254 : const Moose::StateArg * state_limiter) const
255 : {
256 356816 : if (!fi)
257 356816 : fi = _face_info;
258 :
259 356816 : return makeFace(*fi, limiter_type, true, correct_skewness, state_limiter);
260 : }
261 :
262 : bool
263 18633448 : FVFluxKernel::avoidBoundary(const FaceInfo & fi) const
264 : {
265 20501245 : for (const auto bnd_id : fi.boundaryIDs())
266 1868237 : if (_boundaries_to_avoid.count(bnd_id))
267 440 : return true;
268 18633008 : return false;
269 : }
270 :
271 : void
272 0 : FVFluxKernel::computeResidual()
273 : {
274 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
275 : }
276 :
277 : void
278 0 : FVFluxKernel::computeJacobian()
279 : {
280 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
281 : }
282 :
283 : void
284 0 : FVFluxKernel::computeResidualAndJacobian()
285 : {
286 0 : mooseError("FVFluxKernel residual/Jacobian evaluation requires a face information object");
287 : }
288 :
289 : bool
290 12707736 : FVFluxKernel::hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const
291 : {
292 12707736 : if (fi_elem_side)
293 6353868 : return hasBlocks(fi.elem().subdomain_id());
294 : else
295 6353868 : return fi.neighborPtr() && hasBlocks(fi.neighbor().subdomain_id());
296 : }
|
