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 "ArrayDGLowerDKernel.h"
11 :
12 : #include "Assembly.h"
13 : #include "MooseVariable.h"
14 : #include "Problem.h"
15 : #include "SubProblem.h"
16 : #include "SystemBase.h"
17 : #include "MaterialData.h"
18 : #include "ParallelUniqueId.h"
19 :
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/dense_vector.h"
22 : #include "libmesh/numeric_vector.h"
23 : #include "libmesh/dense_subvector.h"
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/quadrature.h"
26 :
27 : InputParameters
28 28787 : ArrayDGLowerDKernel::validParams()
29 : {
30 28787 : InputParameters params = ArrayDGKernel::validParams();
31 28787 : params.addRequiredCoupledVar("lowerd_variable", "Lagrange multiplier");
32 28787 : return params;
33 0 : }
34 :
35 135 : ArrayDGLowerDKernel::ArrayDGLowerDKernel(const InputParameters & parameters)
36 : : ArrayDGKernel(parameters),
37 135 : _lowerd_var(*getArrayVar("lowerd_variable", 0)),
38 135 : _lambda(_is_implicit ? _lowerd_var.slnLower() : _lowerd_var.slnLowerOld()),
39 :
40 135 : _phi_lambda(_lowerd_var.phiLower()),
41 135 : _test_lambda(_lowerd_var.phiLower()),
42 270 : _work_vector(_count)
43 : {
44 135 : const auto & lower_domains = _lowerd_var.activeSubdomains();
45 270 : for (const auto & id : _mesh.interiorLowerDBlocks())
46 135 : if (lower_domains.count(id) == 0)
47 0 : mooseDocumentedError(
48 : "moose",
49 : 29151,
50 0 : "'lowerd_variable' must be defined on the interior lower-dimensional subdomain '" +
51 0 : _mesh.getSubdomainName(id) +
52 : "' that is added by Mesh/build_all_side_lowerd_mesh=true.\nThe check could be overly "
53 : "restrictive.");
54 :
55 270 : for (const auto & id : _var.activeSubdomains())
56 135 : if (_mesh.interiorLowerDBlocks().count(id) > 0)
57 0 : paramError("variable",
58 0 : "Must not be defined on the interior lower-dimensional subdomain '" +
59 0 : _mesh.getSubdomainName(id) + "'");
60 :
61 135 : if (_lowerd_var.count() != _count)
62 0 : paramError("lowerd_variable",
63 : "The number of components must be equal to the number of "
64 : "components of 'variable'");
65 135 : }
66 :
67 : void
68 5364 : ArrayDGLowerDKernel::computeResidual()
69 : {
70 5364 : if (!excludeBoundary())
71 : {
72 5364 : precalculateResidual();
73 :
74 : // Compute the residual for this element
75 5364 : computeElemNeighResidual(Moose::Element);
76 :
77 : // Compute the residual for the neighbor
78 5364 : computeElemNeighResidual(Moose::Neighbor);
79 :
80 5364 : computeLowerDResidual();
81 : }
82 5364 : }
83 :
84 : void
85 5364 : ArrayDGLowerDKernel::computeLowerDResidual()
86 : {
87 5364 : prepareVectorTagLower(_assembly, _lowerd_var.number());
88 :
89 27012 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
90 : {
91 21648 : initLowerDQpResidual();
92 44256 : for (_i = 0; _i < _test_lambda.size(); _i++)
93 : {
94 22608 : _work_vector.setZero();
95 22608 : computeLowerDQpResidual(_work_vector);
96 : mooseAssert(_work_vector.size() == _count,
97 : "Size of local residual is not equal to the number of array variable compoments");
98 22608 : _work_vector *= _JxW[_qp] * _coord[_qp];
99 22608 : _assembly.saveLocalArrayResidual(_local_re, _i, _test_lambda.size(), _work_vector);
100 : }
101 : }
102 :
103 5364 : accumulateTaggedLocalResidual();
104 5364 : }
105 :
106 : void
107 192 : ArrayDGLowerDKernel::computeJacobian()
108 : {
109 192 : if (!excludeBoundary())
110 : {
111 192 : precalculateJacobian();
112 :
113 : // Compute element-element Jacobian
114 192 : computeElemNeighJacobian(Moose::ElementElement);
115 :
116 : // Compute element-neighbor Jacobian
117 192 : computeElemNeighJacobian(Moose::ElementNeighbor);
118 :
119 : // Compute neighbor-element Jacobian
120 192 : computeElemNeighJacobian(Moose::NeighborElement);
121 :
122 : // Compute neighbor-neighbor Jacobian
123 192 : computeElemNeighJacobian(Moose::NeighborNeighbor);
124 :
125 : // Compute the other five pieces of Jacobian related with lower-d variable
126 192 : computeLowerDJacobian(Moose::LowerLower);
127 : }
128 192 : }
129 :
130 : void
131 192 : ArrayDGLowerDKernel::computeLowerDJacobian(Moose::ConstraintJacobianType type)
132 : {
133 : mooseAssert(type != Moose::PrimaryPrimary && type != Moose::PrimarySecondary &&
134 : type != Moose::SecondaryPrimary && type != Moose::SecondarySecondary,
135 : "Jacobian types without lower should be handled in computeElemNeighJacobian");
136 :
137 192 : const ArrayVariableTestValue & test_space =
138 0 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
139 192 : ? _test_lambda
140 : : (type == Moose::PrimaryLower ? _test : _test_neighbor);
141 : unsigned int ivar =
142 0 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
143 192 : ? _lowerd_var.number()
144 0 : : _var.number();
145 :
146 192 : const ArrayVariableTestValue & loc_phi =
147 0 : (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
148 : ? _phi_lambda
149 : : (type == Moose::LowerPrimary ? _phi : _phi_neighbor);
150 : unsigned int jvar =
151 0 : (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
152 192 : ? _lowerd_var.number()
153 0 : : _var.number();
154 :
155 192 : prepareMatrixTagLower(_assembly, ivar, jvar, type);
156 :
157 192 : if (_local_ke.n() == 0 || _local_ke.m() == 0)
158 0 : return;
159 :
160 960 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
161 : {
162 768 : initLowerDQpJacobian(type);
163 1536 : for (_i = 0; _i < test_space.size(); _i++)
164 1536 : for (_j = 0; _j < loc_phi.size(); _j++)
165 : {
166 : // vector or matrix?
167 768 : RealEigenVector v = _JxW[_qp] * _coord[_qp] * computeLowerDQpJacobian(type);
168 768 : _assembly.saveDiagLocalArrayJacobian(
169 768 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), ivar, v);
170 768 : }
171 : }
172 :
173 192 : accumulateTaggedLocalMatrix();
174 : }
175 :
176 : void
177 1920 : ArrayDGLowerDKernel::computeOffDiagJacobian(const unsigned int jvar_num)
178 : {
179 1920 : if (!excludeBoundary())
180 : {
181 1920 : precalculateOffDiagJacobian(jvar_num);
182 :
183 : /*
184 : * Note: we cannot call compute jacobian functions like in DGLowerDKernel
185 : * because we could have cross component couplings for the array variables
186 : */
187 :
188 1920 : const auto & jvar = getVariable(jvar_num);
189 :
190 : // Compute element-element Jacobian
191 1920 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
192 :
193 : // Compute element-neighbor Jacobian
194 1920 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
195 :
196 : // Compute neighbor-element Jacobian
197 1920 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
198 :
199 : // Compute neighbor-neighbor Jacobian
200 1920 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
201 :
202 : // Compute the other five pieces of Jacobian related with lower-d variable
203 1920 : computeOffDiagLowerDJacobian(Moose::LowerLower, jvar);
204 1920 : computeOffDiagLowerDJacobian(Moose::LowerSecondary, jvar);
205 1920 : computeOffDiagLowerDJacobian(Moose::LowerPrimary, jvar);
206 1920 : computeOffDiagLowerDJacobian(Moose::SecondaryLower, jvar);
207 1920 : computeOffDiagLowerDJacobian(Moose::PrimaryLower, jvar);
208 : }
209 1920 : }
210 :
211 : void
212 9600 : ArrayDGLowerDKernel::computeOffDiagLowerDJacobian(Moose::ConstraintJacobianType type,
213 : const MooseVariableFEBase & jvar)
214 : {
215 : mooseAssert(type != Moose::PrimaryPrimary && type != Moose::PrimarySecondary &&
216 : type != Moose::SecondaryPrimary && type != Moose::SecondarySecondary,
217 : "Jacobian types without lower should be handled in computeOffDiagElemNeighJacobian");
218 :
219 9600 : const ArrayVariableTestValue & test_space =
220 7680 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
221 17280 : ? _test_lambda
222 : : (type == Moose::PrimaryLower ? _test : _test_neighbor);
223 : unsigned int ivar =
224 7680 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
225 13440 : ? _lowerd_var.number()
226 3840 : : _var.number();
227 :
228 9600 : prepareMatrixTagLower(_assembly, ivar, jvar.number(), type);
229 9600 : if (_local_ke.n() == 0 || _local_ke.m() == 0)
230 4800 : return;
231 :
232 4800 : if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
233 : {
234 0 : const auto & jv0 = static_cast<const MooseVariable &>(jvar);
235 : const VariableTestValue & loc_phi =
236 0 : (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
237 0 : ? jv0.phiLower()
238 0 : : (type == Moose::LowerPrimary ? jv0.phiFace() : jv0.phiFaceNeighbor());
239 :
240 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
241 : {
242 0 : initLowerDQpOffDiagJacobian(type, jvar);
243 0 : for (_i = 0; _i < test_space.size(); _i++)
244 0 : for (_j = 0; _j < loc_phi.size(); _j++)
245 : {
246 0 : RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeLowerDQpOffDiagJacobian(type, jvar);
247 0 : _assembly.saveFullLocalArrayJacobian(
248 0 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), ivar, jvar.number(), v);
249 0 : }
250 : }
251 : }
252 4800 : else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
253 : {
254 4800 : const auto & jv1 = static_cast<const ArrayMooseVariable &>(jvar);
255 : const ArrayVariableTestValue & loc_phi =
256 3840 : (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
257 2880 : ? jv1.phiLower()
258 6720 : : (type == Moose::LowerPrimary ? jv1.phiFace() : jv1.phiFaceNeighbor());
259 :
260 24480 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
261 : {
262 19680 : initLowerDQpOffDiagJacobian(type, jvar);
263 116448 : for (_i = 0; _i < test_space.size(); _i++)
264 297984 : for (_j = 0; _j < loc_phi.size(); _j++)
265 : {
266 201216 : RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeLowerDQpOffDiagJacobian(type, jvar);
267 201216 : _assembly.saveFullLocalArrayJacobian(
268 201216 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), ivar, jvar.number(), v);
269 201216 : }
270 : }
271 : }
272 : else
273 0 : mooseError("Vector variable cannot be coupled into array DG kernel currently");
274 :
275 4800 : accumulateTaggedLocalMatrix();
276 : }
277 :
278 : RealEigenMatrix
279 185472 : ArrayDGLowerDKernel::computeLowerDQpOffDiagJacobian(Moose::ConstraintJacobianType type,
280 : const MooseVariableFEBase & jvar)
281 : {
282 185472 : if (jvar.number() == _var.number())
283 : {
284 90240 : if (type == Moose::LowerSecondary || type == Moose::LowerPrimary)
285 180480 : return computeLowerDQpJacobian(type).asDiagonal();
286 : }
287 95232 : else if (jvar.number() == _lowerd_var.number())
288 : {
289 95232 : if (type == Moose::SecondaryLower || type == Moose::PrimaryLower || type == Moose::LowerLower)
290 190464 : return computeLowerDQpJacobian(type).asDiagonal();
291 : }
292 :
293 0 : return RealEigenMatrix::Zero(_var.count(), jvar.count());
294 : }
|
