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 "DGKernelBase.h"
11 : #include "MooseVariable.h"
12 : #include "Problem.h"
13 : #include "SubProblem.h"
14 : #include "SystemBase.h"
15 : #include "MaterialData.h"
16 : #include "ParallelUniqueId.h"
17 :
18 : #include "libmesh/dof_map.h"
19 : #include "libmesh/dense_vector.h"
20 : #include "libmesh/numeric_vector.h"
21 : #include "libmesh/dense_subvector.h"
22 : #include "libmesh/libmesh_common.h"
23 : #include "libmesh/quadrature.h"
24 :
25 : InputParameters
26 45557 : DGKernelBase::validParams()
27 : {
28 45557 : InputParameters params = NeighborResidualObject::validParams();
29 45557 : params += TwoMaterialPropertyInterface::validParams();
30 45557 : params += BlockRestrictable::validParams();
31 45557 : params += BoundaryRestrictable::validParams();
32 136671 : params.addParam<bool>("use_displaced_mesh",
33 91114 : false,
34 : "Whether or not this object should use the "
35 : "displaced mesh for computation. Note that in "
36 : "the case this is true but no displacements "
37 : "are provided in the Mesh block the "
38 : "undisplaced mesh will still be used.");
39 182228 : params.addParamNamesToGroup("use_displaced_mesh", "Advanced");
40 :
41 182228 : params.addParam<std::vector<AuxVariableName>>(
42 : "save_in",
43 : {},
44 : "The name of auxiliary variables to save this Kernel's residual contributions to. "
45 : " Everything about that variable must match everything about this variable (the "
46 : "type, what blocks it's on, etc.)");
47 182228 : params.addParam<std::vector<AuxVariableName>>(
48 : "diag_save_in",
49 : {},
50 : "The name of auxiliary variables to save this Kernel's diagonal Jacobian "
51 : "contributions to. Everything about that variable must match everything "
52 : "about this variable (the type, what blocks it's on, etc.)");
53 136671 : params.addParamNamesToGroup("diag_save_in save_in", "Advanced");
54 :
55 : // DG Kernels always need one layer of ghosting.
56 136671 : params.addRelationshipManager("ElementSideNeighborLayers",
57 : Moose::RelationshipManagerType::GEOMETRIC |
58 : Moose::RelationshipManagerType::ALGEBRAIC |
59 : Moose::RelationshipManagerType::COUPLING);
60 :
61 182228 : params.addParam<std::vector<BoundaryName>>(
62 : "exclude_boundary", "The internal side sets to be excluded from this kernel.");
63 91114 : params.registerBase("DGKernel");
64 45557 : params.registerSystemAttributeName("DGKernel");
65 :
66 45557 : return params;
67 0 : }
68 :
69 1419 : DGKernelBase::DGKernelBase(const InputParameters & parameters)
70 : : NeighborResidualObject(parameters),
71 : BlockRestrictable(this),
72 : BoundaryRestrictable(this, false), // false for _not_ nodal
73 : NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, false, false),
74 : TwoMaterialPropertyInterface(this, blockIDs(), boundaryIDs()),
75 : ElementIDInterface(this),
76 2838 : _current_elem(_assembly.elem()),
77 1419 : _current_elem_volume(_assembly.elemVolume()),
78 :
79 1419 : _neighbor_elem(_assembly.neighbor()),
80 1419 : _neighbor_elem_volume(_assembly.neighborVolume()),
81 :
82 1419 : _current_side(_assembly.side()),
83 1419 : _current_side_elem(_assembly.sideElem()),
84 1419 : _current_side_volume(_assembly.sideElemVolume()),
85 :
86 1419 : _coord_sys(_assembly.coordSystem()),
87 1419 : _q_point(_assembly.qPointsFace()),
88 1419 : _qrule(_assembly.qRuleFace()),
89 1419 : _JxW(_assembly.JxWFace()),
90 1419 : _coord(_assembly.coordTransformation()),
91 1419 : _normals(_assembly.normals()),
92 :
93 1419 : _save_in_strings(parameters.get<std::vector<AuxVariableName>>("save_in")),
94 4257 : _diag_save_in_strings(parameters.get<std::vector<AuxVariableName>>("diag_save_in"))
95 : {
96 : // Gather information on broken boundaries
97 2838 : std::vector<BoundaryName> bnd = isParamValid("exclude_boundary")
98 1419 : ? getParam<std::vector<BoundaryName>>("exclude_boundary")
99 4257 : : std::vector<BoundaryName>(0);
100 1419 : auto bnd_ids = _mesh.getBoundaryIDs(bnd);
101 :
102 : // check if the broken boundary ids are valid
103 1419 : auto & valid_ids = _mesh.meshSidesetIds();
104 1419 : std::vector<BoundaryName> diff;
105 1419 : for (unsigned int i = 0; i < bnd_ids.size(); ++i)
106 0 : if (valid_ids.find(bnd_ids[i]) == valid_ids.end())
107 0 : diff.push_back(bnd[i]);
108 1419 : if (!diff.empty())
109 : {
110 0 : auto msg = "DGKernel '" + name() +
111 : "' contains the following boundary names that do not exist on the mesh: " +
112 0 : Moose::stringify(diff, ", ");
113 0 : paramError("exclude_boundary", msg);
114 0 : }
115 :
116 1419 : _excluded_boundaries.insert(bnd_ids.begin(), bnd_ids.end());
117 1419 : }
118 :
119 : void
120 2588071 : DGKernelBase::computeResidual()
121 : {
122 2588071 : if (!excludeBoundary())
123 : {
124 2588071 : precalculateResidual();
125 :
126 : // Compute the residual for this element
127 2588071 : computeElemNeighResidual(Moose::Element);
128 :
129 : // Compute the residual for the neighbor
130 2588071 : computeElemNeighResidual(Moose::Neighbor);
131 : }
132 2588071 : }
133 :
134 : void
135 99626 : DGKernelBase::computeJacobian()
136 : {
137 99626 : if (!excludeBoundary())
138 : {
139 99626 : precalculateJacobian();
140 :
141 : // Compute element-element Jacobian
142 99626 : computeElemNeighJacobian(Moose::ElementElement);
143 :
144 : // Compute element-neighbor Jacobian
145 99626 : computeElemNeighJacobian(Moose::ElementNeighbor);
146 :
147 : // Compute neighbor-element Jacobian
148 99626 : computeElemNeighJacobian(Moose::NeighborElement);
149 :
150 : // Compute neighbor-neighbor Jacobian
151 99626 : computeElemNeighJacobian(Moose::NeighborNeighbor);
152 : }
153 99626 : }
154 :
155 : void
156 18168 : DGKernelBase::computeOffDiagJacobian(const unsigned int jvar_num)
157 : {
158 18168 : if (!excludeBoundary())
159 : {
160 18168 : const auto & jvar = getVariable(jvar_num);
161 :
162 18168 : if (jvar_num == variable().number())
163 15210 : computeJacobian();
164 : else
165 : {
166 2958 : precalculateOffDiagJacobian(jvar_num);
167 :
168 : // Compute element-element Jacobian
169 2958 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
170 :
171 : // Compute element-neighbor Jacobian
172 2958 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
173 :
174 : // Compute neighbor-element Jacobian
175 2958 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
176 :
177 : // Compute neighbor-neighbor Jacobian
178 2958 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
179 : }
180 : }
181 18168 : }
182 :
183 : bool
184 2760597 : DGKernelBase::excludeBoundary() const
185 : {
186 2760597 : if (_excluded_boundaries.empty())
187 2760597 : return false;
188 :
189 0 : auto boundary_ids = _mesh.getBoundaryIDs(_current_elem, _current_side);
190 0 : for (auto bid : boundary_ids)
191 0 : if (_excluded_boundaries.find(bid) != _excluded_boundaries.end())
192 0 : return true;
193 :
194 : // make sure we will also break on the neighboring side
195 0 : unsigned int neighbor_side = _neighbor_elem->which_neighbor_am_i(_current_elem);
196 0 : boundary_ids = _mesh.getBoundaryIDs(_neighbor_elem, neighbor_side);
197 0 : for (auto bid : boundary_ids)
198 0 : if (_excluded_boundaries.find(bid) != _excluded_boundaries.end())
199 0 : return true;
200 :
201 0 : return false;
202 0 : }
203 :
204 : void
205 149416 : DGKernelBase::prepareShapes(const unsigned int var_num)
206 : {
207 149416 : _subproblem.prepareFaceShapes(var_num, _tid);
208 149416 : }
|
