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 28807 : ArrayDGLowerDKernel::validParams()
29 : {
30 28807 : InputParameters params = ArrayDGKernel::validParams();
31 28807 : params.addRequiredCoupledVar("lowerd_variable", "Lagrange multiplier");
32 28807 : return params;
33 0 : }
34 :
35 145 : ArrayDGLowerDKernel::ArrayDGLowerDKernel(const InputParameters & parameters)
36 : : ArrayDGKernel(parameters),
37 145 : _lowerd_var(*getArrayVar("lowerd_variable", 0)),
38 145 : _lambda(_is_implicit ? _lowerd_var.slnLower() : _lowerd_var.slnLowerOld()),
39 :
40 145 : _phi_lambda(_lowerd_var.phiLower()),
41 145 : _test_lambda(_lowerd_var.phiLower()),
42 290 : _work_vector(_count)
43 : {
44 145 : const auto & lower_domains = _lowerd_var.activeSubdomains();
45 290 : for (const auto & id : _mesh.interiorLowerDBlocks())
46 145 : 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 290 : for (const auto & id : _var.activeSubdomains())
56 145 : 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 145 : 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 145 : }
66 :
67 : void
68 6036 : ArrayDGLowerDKernel::computeResidual()
69 : {
70 6036 : if (!excludeBoundary())
71 : {
72 6036 : precalculateResidual();
73 :
74 : // Compute the residual for this element
75 6036 : computeElemNeighResidual(Moose::Element);
76 :
77 : // Compute the residual for the neighbor
78 6036 : computeElemNeighResidual(Moose::Neighbor);
79 :
80 6036 : computeLowerDResidual();
81 : }
82 6036 : }
83 :
84 : void
85 6036 : ArrayDGLowerDKernel::computeLowerDResidual()
86 : {
87 6036 : prepareVectorTagLower(_assembly, _lowerd_var.number());
88 :
89 30396 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
90 : {
91 24360 : initLowerDQpResidual();
92 49800 : for (_i = 0; _i < _test_lambda.size(); _i++)
93 : {
94 25440 : _work_vector.setZero();
95 25440 : 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 25440 : _work_vector *= _JxW[_qp] * _coord[_qp];
99 25440 : _assembly.saveLocalArrayResidual(_local_re, _i, _test_lambda.size(), _work_vector);
100 : }
101 : }
102 :
103 6036 : accumulateTaggedLocalResidual();
104 6036 : }
105 :
106 : void
107 216 : ArrayDGLowerDKernel::computeJacobian()
108 : {
109 216 : if (!excludeBoundary())
110 : {
111 216 : precalculateJacobian();
112 :
113 : // Compute element-element Jacobian
114 216 : computeElemNeighJacobian(Moose::ElementElement);
115 :
116 : // Compute element-neighbor Jacobian
117 216 : computeElemNeighJacobian(Moose::ElementNeighbor);
118 :
119 : // Compute neighbor-element Jacobian
120 216 : computeElemNeighJacobian(Moose::NeighborElement);
121 :
122 : // Compute neighbor-neighbor Jacobian
123 216 : computeElemNeighJacobian(Moose::NeighborNeighbor);
124 :
125 : // Compute the other five pieces of Jacobian related with lower-d variable
126 216 : computeLowerDJacobian(Moose::LowerLower);
127 : }
128 216 : }
129 :
130 : void
131 216 : 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 216 : const ArrayVariableTestValue & test_space =
138 0 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
139 216 : ? _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 216 : ? _lowerd_var.number()
144 0 : : _var.number();
145 :
146 216 : 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 216 : ? _lowerd_var.number()
153 0 : : _var.number();
154 :
155 216 : prepareMatrixTagLower(_assembly, ivar, jvar, type);
156 :
157 216 : if (_local_ke.n() == 0 || _local_ke.m() == 0)
158 0 : return;
159 :
160 1080 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
161 : {
162 864 : initLowerDQpJacobian(type);
163 1728 : for (_i = 0; _i < test_space.size(); _i++)
164 1728 : for (_j = 0; _j < loc_phi.size(); _j++)
165 : {
166 : // vector or matrix?
167 864 : RealEigenVector v = _JxW[_qp] * _coord[_qp] * computeLowerDQpJacobian(type);
168 864 : _assembly.saveDiagLocalArrayJacobian(
169 864 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), ivar, v);
170 864 : }
171 : }
172 :
173 216 : accumulateTaggedLocalMatrix();
174 : }
175 :
176 : void
177 2160 : ArrayDGLowerDKernel::computeOffDiagJacobian(const unsigned int jvar_num)
178 : {
179 2160 : if (!excludeBoundary())
180 : {
181 2160 : 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 2160 : const auto & jvar = getVariable(jvar_num);
189 :
190 : // Compute element-element Jacobian
191 2160 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
192 :
193 : // Compute element-neighbor Jacobian
194 2160 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
195 :
196 : // Compute neighbor-element Jacobian
197 2160 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
198 :
199 : // Compute neighbor-neighbor Jacobian
200 2160 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
201 :
202 : // Compute the other five pieces of Jacobian related with lower-d variable
203 2160 : computeOffDiagLowerDJacobian(Moose::LowerLower, jvar);
204 2160 : computeOffDiagLowerDJacobian(Moose::LowerSecondary, jvar);
205 2160 : computeOffDiagLowerDJacobian(Moose::LowerPrimary, jvar);
206 2160 : computeOffDiagLowerDJacobian(Moose::SecondaryLower, jvar);
207 2160 : computeOffDiagLowerDJacobian(Moose::PrimaryLower, jvar);
208 : }
209 2160 : }
210 :
211 : void
212 10800 : 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 10800 : const ArrayVariableTestValue & test_space =
220 8640 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
221 19440 : ? _test_lambda
222 : : (type == Moose::PrimaryLower ? _test : _test_neighbor);
223 : unsigned int ivar =
224 8640 : (type == Moose::LowerLower || type == Moose::LowerSecondary || type == Moose::LowerPrimary)
225 15120 : ? _lowerd_var.number()
226 4320 : : _var.number();
227 :
228 10800 : prepareMatrixTagLower(_assembly, ivar, jvar.number(), type);
229 10800 : if (_local_ke.n() == 0 || _local_ke.m() == 0)
230 5400 : return;
231 :
232 5400 : 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 5400 : else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
253 : {
254 5400 : const auto & jv1 = static_cast<const ArrayMooseVariable &>(jvar);
255 : const ArrayVariableTestValue & loc_phi =
256 4320 : (type == Moose::LowerLower || type == Moose::SecondaryLower || type == Moose::PrimaryLower)
257 3240 : ? jv1.phiLower()
258 7560 : : (type == Moose::LowerPrimary ? jv1.phiFace() : jv1.phiFaceNeighbor());
259 :
260 27540 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
261 : {
262 22140 : initLowerDQpOffDiagJacobian(type, jvar);
263 131004 : for (_i = 0; _i < test_space.size(); _i++)
264 335232 : for (_j = 0; _j < loc_phi.size(); _j++)
265 : {
266 226368 : RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeLowerDQpOffDiagJacobian(type, jvar);
267 226368 : _assembly.saveFullLocalArrayJacobian(
268 226368 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), ivar, jvar.number(), v);
269 226368 : }
270 : }
271 : }
272 : else
273 0 : mooseError("Vector variable cannot be coupled into array DG kernel currently");
274 :
275 5400 : accumulateTaggedLocalMatrix();
276 : }
277 :
278 : RealEigenMatrix
279 208656 : ArrayDGLowerDKernel::computeLowerDQpOffDiagJacobian(Moose::ConstraintJacobianType type,
280 : const MooseVariableFEBase & jvar)
281 : {
282 208656 : if (jvar.number() == _var.number())
283 : {
284 101520 : if (type == Moose::LowerSecondary || type == Moose::LowerPrimary)
285 203040 : return computeLowerDQpJacobian(type).asDiagonal();
286 : }
287 107136 : else if (jvar.number() == _lowerd_var.number())
288 : {
289 107136 : if (type == Moose::SecondaryLower || type == Moose::PrimaryLower || type == Moose::LowerLower)
290 214272 : return computeLowerDQpJacobian(type).asDiagonal();
291 : }
292 :
293 0 : return RealEigenMatrix::Zero(_var.count(), jvar.count());
294 : }
|
