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 "ArrayDGKernel.h"
11 : #include "Assembly.h"
12 : #include "MooseVariable.h"
13 : #include "Problem.h"
14 : #include "SubProblem.h"
15 : #include "SystemBase.h"
16 : #include "MaterialData.h"
17 : #include "ParallelUniqueId.h"
18 :
19 : #include "libmesh/dof_map.h"
20 : #include "libmesh/dense_vector.h"
21 : #include "libmesh/numeric_vector.h"
22 : #include "libmesh/dense_subvector.h"
23 : #include "libmesh/libmesh_common.h"
24 : #include "libmesh/quadrature.h"
25 :
26 : InputParameters
27 9520 : ArrayDGKernel::validParams()
28 : {
29 9520 : InputParameters params = DGKernelBase::validParams();
30 9520 : return params;
31 : }
32 :
33 178 : ArrayDGKernel::ArrayDGKernel(const InputParameters & parameters)
34 : : DGKernelBase(parameters),
35 : NeighborMooseVariableInterface(
36 : this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_ARRAY),
37 356 : _var(*mooseVariable()),
38 178 : _u(_is_implicit ? _var.sln() : _var.slnOld()),
39 178 : _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
40 :
41 178 : _phi(_var.phiFace()),
42 178 : _grad_phi(_var.gradPhiFace()),
43 178 : _array_grad_phi(_var.arrayGradPhiFace()),
44 :
45 178 : _test(_var.phiFace()),
46 178 : _grad_test(_var.gradPhiFace()),
47 178 : _array_grad_test(_var.arrayGradPhiFace()),
48 :
49 178 : _phi_neighbor(_var.phiFaceNeighbor()),
50 178 : _grad_phi_neighbor(_var.gradPhiFaceNeighbor()),
51 178 : _array_grad_phi_neighbor(_var.arrayGradPhiFaceNeighbor()),
52 :
53 178 : _test_neighbor(_var.phiFaceNeighbor()),
54 178 : _grad_test_neighbor(_var.gradPhiFaceNeighbor()),
55 178 : _array_grad_test_neighbor(_var.arrayGradPhiFaceNeighbor()),
56 :
57 178 : _u_neighbor(_is_implicit ? _var.slnNeighbor() : _var.slnOldNeighbor()),
58 178 : _grad_u_neighbor(_is_implicit ? _var.gradSlnNeighbor() : _var.gradSlnOldNeighbor()),
59 :
60 178 : _array_normals(_assembly.mappedNormals()),
61 178 : _count(_var.count()),
62 :
63 356 : _work_vector(_count)
64 : {
65 178 : addMooseVariableDependency(mooseVariable());
66 :
67 178 : _save_in.resize(_save_in_strings.size());
68 178 : _diag_save_in.resize(_diag_save_in_strings.size());
69 :
70 191 : for (unsigned int i = 0; i < _save_in_strings.size(); i++)
71 : {
72 13 : MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
73 13 : _save_in_strings[i],
74 : Moose::VarKindType::VAR_AUXILIARY,
75 : Moose::VarFieldType::VAR_FIELD_ARRAY);
76 :
77 13 : if (_sys.hasVariable(_save_in_strings[i]))
78 0 : mooseError("Trying to use solution variable " + _save_in_strings[i] +
79 0 : " as a save_in variable in " + name());
80 :
81 13 : if (var->feType() != _var.feType())
82 0 : paramError(
83 : "save_in",
84 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
85 0 : moose::internal::incompatVarMsg(*var, _var));
86 :
87 13 : _save_in[i] = var;
88 13 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
89 13 : addMooseVariableDependency(var);
90 : }
91 :
92 178 : _has_save_in = _save_in.size() > 0;
93 :
94 191 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
95 : {
96 13 : MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
97 13 : _diag_save_in_strings[i],
98 : Moose::VarKindType::VAR_AUXILIARY,
99 : Moose::VarFieldType::VAR_FIELD_ARRAY);
100 :
101 13 : if (_sys.hasVariable(_diag_save_in_strings[i]))
102 0 : mooseError("Trying to use solution variable " + _diag_save_in_strings[i] +
103 0 : " as a diag_save_in variable in " + name());
104 :
105 13 : if (var->feType() != _var.feType())
106 0 : paramError(
107 : "diag_save_in",
108 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
109 0 : moose::internal::incompatVarMsg(*var, _var));
110 :
111 13 : _diag_save_in[i] = var;
112 13 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
113 13 : addMooseVariableDependency(var);
114 : }
115 :
116 178 : _has_diag_save_in = _diag_save_in.size() > 0;
117 178 : }
118 :
119 : void
120 16680 : ArrayDGKernel::computeElemNeighResidual(Moose::DGResidualType type)
121 : {
122 : bool is_elem;
123 16680 : if (type == Moose::Element)
124 8340 : is_elem = true;
125 : else
126 8340 : is_elem = false;
127 :
128 16680 : const ArrayVariableTestValue & test_space = is_elem ? _test : _test_neighbor;
129 :
130 16680 : if (is_elem)
131 8340 : prepareVectorTag(_assembly, _var.number());
132 : else
133 8340 : prepareVectorTagNeighbor(_assembly, _var.number());
134 :
135 71880 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
136 : {
137 55200 : initQpResidual(type);
138 545376 : for (_i = 0; _i < test_space.size(); _i++)
139 : {
140 490176 : _work_vector.setZero();
141 490176 : computeQpResidual(type, _work_vector);
142 : mooseAssert(_work_vector.size() == _count,
143 : "Size of local residual is not equal to the number of array variable compoments");
144 490176 : _work_vector *= _JxW[_qp] * _coord[_qp];
145 490176 : _assembly.saveLocalArrayResidual(_local_re, _i, test_space.size(), _work_vector);
146 : }
147 : }
148 :
149 16680 : accumulateTaggedLocalResidual();
150 :
151 16680 : if (_has_save_in)
152 2304 : for (const auto & var : _save_in)
153 : {
154 1152 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
155 1152 : if (!avar)
156 0 : mooseError("Save-in variable for an array kernel must be an array variable");
157 :
158 1152 : if (is_elem)
159 576 : avar->addSolution(_local_re);
160 : else
161 576 : avar->addSolutionNeighbor(_local_re);
162 : }
163 16680 : }
164 :
165 : void
166 768 : ArrayDGKernel::computeElemNeighJacobian(Moose::DGJacobianType type)
167 : {
168 768 : const ArrayVariableTestValue & test_space =
169 768 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
170 1344 : const ArrayVariableTestValue & loc_phi =
171 576 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? _phi : _phi_neighbor;
172 :
173 768 : if (type == Moose::ElementElement)
174 192 : prepareMatrixTag(_assembly, _var.number(), _var.number());
175 : else
176 576 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), type);
177 :
178 3840 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
179 : {
180 3072 : initQpJacobian(type);
181 33792 : for (_i = 0; _i < test_space.size(); _i++)
182 337920 : for (_j = 0; _j < loc_phi.size(); _j++)
183 : {
184 307200 : RealEigenVector v = _JxW[_qp] * _coord[_qp] * computeQpJacobian(type);
185 307200 : _assembly.saveDiagLocalArrayJacobian(
186 614400 : _local_ke, _i, test_space.size(), _j, loc_phi.size(), _var.number(), v);
187 307200 : }
188 : }
189 :
190 768 : accumulateTaggedLocalMatrix();
191 :
192 768 : if (_has_diag_save_in && (type == Moose::ElementElement || type == Moose::NeighborNeighbor))
193 : {
194 0 : DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
195 0 : for (const auto & var : _diag_save_in)
196 : {
197 0 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
198 0 : if (!avar)
199 0 : mooseError("Save-in variable for an array kernel must be an array variable");
200 :
201 0 : if (type == Moose::ElementElement)
202 0 : avar->addSolution(diag);
203 : else
204 0 : avar->addSolutionNeighbor(diag);
205 : }
206 0 : }
207 768 : }
208 :
209 : RealEigenVector
210 2121600 : ArrayDGKernel::computeQpJacobian(Moose::DGJacobianType)
211 : {
212 2121600 : return RealEigenVector::Zero(_count);
213 : }
214 :
215 : void
216 1344 : ArrayDGKernel::computeOffDiagJacobian(const unsigned int jvar_num)
217 : {
218 1344 : if (!excludeBoundary())
219 : {
220 1344 : const auto & jvar = getVariable(jvar_num);
221 :
222 1344 : precalculateOffDiagJacobian(jvar_num);
223 :
224 : // Compute element-element Jacobian
225 1344 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
226 :
227 : // Compute element-neighbor Jacobian
228 1344 : computeOffDiagElemNeighJacobian(Moose::ElementNeighbor, jvar);
229 :
230 : // Compute neighbor-element Jacobian
231 1344 : computeOffDiagElemNeighJacobian(Moose::NeighborElement, jvar);
232 :
233 : // Compute neighbor-neighbor Jacobian
234 1344 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
235 : }
236 1344 : }
237 :
238 : void
239 13056 : ArrayDGKernel::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type,
240 : const MooseVariableFEBase & jvar)
241 : {
242 13056 : const ArrayVariableTestValue & test_space =
243 13056 : (type == Moose::ElementElement || type == Moose::ElementNeighbor) ? _test : _test_neighbor;
244 :
245 13056 : if (type == Moose::ElementElement)
246 3264 : prepareMatrixTag(_assembly, _var.number(), jvar.number());
247 : else
248 9792 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar.number(), type);
249 :
250 13056 : if (_local_ke.n() == 0 || _local_ke.m() == 0)
251 3840 : return;
252 :
253 9216 : if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD)
254 : {
255 0 : const auto & jv0 = static_cast<const MooseVariable &>(jvar);
256 : const VariableTestValue & loc_phi =
257 0 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? jv0.phiFace()
258 0 : : jv0.phiFaceNeighbor();
259 :
260 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
261 : {
262 0 : initQpOffDiagJacobian(type, jvar);
263 0 : for (_i = 0; _i < test_space.size(); _i++)
264 0 : for (_j = 0; _j < loc_phi.size(); _j++)
265 : {
266 0 : RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
267 0 : _assembly.saveFullLocalArrayJacobian(_local_ke,
268 : _i,
269 : test_space.size(),
270 : _j,
271 : loc_phi.size(),
272 0 : _var.number(),
273 : jvar.number(),
274 : v);
275 0 : }
276 : }
277 : }
278 9216 : else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY)
279 : {
280 9216 : const auto & jv1 = static_cast<const ArrayMooseVariable &>(jvar);
281 : const ArrayVariableTestValue & loc_phi =
282 9216 : (type == Moose::ElementElement || type == Moose::NeighborElement) ? jv1.phiFace()
283 13824 : : jv1.phiFaceNeighbor();
284 :
285 35712 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
286 : {
287 26496 : initQpOffDiagJacobian(type, jvar);
288 236544 : for (_i = 0; _i < test_space.size(); _i++)
289 2196480 : for (_j = 0; _j < loc_phi.size(); _j++)
290 : {
291 1986432 : RealEigenMatrix v = _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(type, jvar);
292 3972864 : _assembly.saveFullLocalArrayJacobian(_local_ke,
293 : _i,
294 : test_space.size(),
295 : _j,
296 : loc_phi.size(),
297 1986432 : _var.number(),
298 : jvar.number(),
299 : v);
300 1986432 : }
301 : }
302 : }
303 : else
304 0 : mooseError("Vector variable cannot be coupled into array DG kernel currently");
305 :
306 9216 : accumulateTaggedLocalMatrix();
307 :
308 9984 : if (_has_diag_save_in && (type == Moose::ElementElement || type == Moose::NeighborNeighbor) &&
309 768 : _var.number() == jvar.number())
310 : {
311 768 : DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
312 1536 : for (const auto & var : _diag_save_in)
313 : {
314 768 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
315 768 : if (!avar)
316 0 : mooseError("Save-in variable for an array kernel must be an array variable");
317 :
318 768 : if (type == Moose::ElementElement)
319 384 : avar->addSolution(diag);
320 : else
321 384 : avar->addSolutionNeighbor(diag);
322 : }
323 768 : }
324 : }
325 :
326 : RealEigenMatrix
327 1986432 : ArrayDGKernel::computeQpOffDiagJacobian(Moose::DGJacobianType type,
328 : const MooseVariableFEBase & jvar)
329 : {
330 1986432 : if (jvar.number() == _var.number())
331 1986432 : return computeQpJacobian(type).asDiagonal();
332 : else
333 0 : return RealEigenMatrix::Zero(_var.count(), jvar.count());
334 : }
|
