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 "ArrayKernel.h"
11 :
12 : #include "Assembly.h"
13 : #include "MooseVariableFE.h"
14 : #include "MooseVariableScalar.h"
15 : #include "SubProblem.h"
16 : #include "NonlinearSystem.h"
17 : #include "FEProblemBase.h"
18 :
19 : #include "libmesh/threads.h"
20 : #include "libmesh/quadrature.h"
21 :
22 : InputParameters
23 103430 : ArrayKernel::validParams()
24 : {
25 103430 : InputParameters params = KernelBase::validParams();
26 103430 : params.registerBase("ArrayKernel");
27 103430 : return params;
28 0 : }
29 :
30 1835 : ArrayKernel::ArrayKernel(const InputParameters & parameters)
31 : : KernelBase(parameters),
32 : MooseVariableInterface<RealEigenVector>(this,
33 : false,
34 : "variable",
35 : Moose::VarKindType::VAR_SOLVER,
36 : Moose::VarFieldType::VAR_FIELD_ARRAY),
37 3670 : _var(*mooseVariable()),
38 1835 : _test(_var.phi()),
39 1835 : _grad_test(_var.gradPhi()),
40 1835 : _array_grad_test(_var.arrayGradPhi()),
41 1835 : _phi(_assembly.phi(_var)),
42 1835 : _grad_phi(_assembly.gradPhi(_var)),
43 1835 : _u(_is_implicit ? _var.sln() : _var.slnOld()),
44 1835 : _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
45 1835 : _count(_var.count()),
46 3670 : _work_vector(_count)
47 : {
48 1835 : addMooseVariableDependency(mooseVariable());
49 :
50 1835 : _save_in.resize(_save_in_strings.size());
51 1835 : _diag_save_in.resize(_diag_save_in_strings.size());
52 :
53 1848 : for (unsigned int i = 0; i < _save_in_strings.size(); i++)
54 : {
55 13 : ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _save_in_strings[i]);
56 :
57 13 : if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_save_in_strings[i]))
58 0 : paramError("save_in", "cannot use solution variable as save-in variable");
59 :
60 13 : if (var->feType() != _var.feType())
61 0 : paramError(
62 : "save_in",
63 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
64 0 : moose::internal::incompatVarMsg(*var, _var));
65 :
66 13 : _save_in[i] = var;
67 13 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
68 13 : addMooseVariableDependency(var);
69 : }
70 :
71 1835 : _has_save_in = _save_in.size() > 0;
72 :
73 1848 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
74 : {
75 13 : ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _diag_save_in_strings[i]);
76 :
77 13 : if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_diag_save_in_strings[i]))
78 0 : paramError("diag_save_in", "cannot use solution variable as diag save-in variable");
79 :
80 13 : if (var->feType() != _var.feType())
81 0 : paramError(
82 : "diag_save_in",
83 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
84 0 : moose::internal::incompatVarMsg(*var, _var));
85 :
86 13 : _diag_save_in[i] = var;
87 13 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
88 13 : addMooseVariableDependency(var);
89 : }
90 :
91 1835 : _has_diag_save_in = _diag_save_in.size() > 0;
92 1835 : }
93 :
94 : void
95 4595757 : ArrayKernel::computeResidual()
96 : {
97 4595757 : prepareVectorTag(_assembly, _var.number());
98 :
99 4595757 : precalculateResidual();
100 26045877 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
101 : {
102 21450120 : initQpResidual();
103 131555256 : for (_i = 0; _i < _test.size(); _i++)
104 : {
105 110105136 : _work_vector.setZero();
106 110105136 : computeQpResidual(_work_vector);
107 : mooseAssert(_work_vector.size() == _count,
108 : "Size of local residual is not equal to the number of array variable compoments");
109 110105136 : _work_vector *= _JxW[_qp] * _coord[_qp];
110 110105136 : _assembly.saveLocalArrayResidual(_local_re, _i, _test.size(), _work_vector);
111 : }
112 : }
113 :
114 4595757 : accumulateTaggedLocalResidual();
115 :
116 4595757 : if (_has_save_in)
117 : {
118 384 : libMesh::Threads::spin_mutex::scoped_lock lock(libMesh::Threads::spin_mtx);
119 768 : for (const auto & var : _save_in)
120 : {
121 384 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
122 384 : if (avar)
123 384 : avar->addSolution(_local_re);
124 : else
125 0 : mooseError("Save-in variable for an array kernel must be an array variable");
126 : }
127 384 : }
128 4595757 : }
129 :
130 : void
131 458648 : ArrayKernel::computeJacobian()
132 : {
133 458648 : prepareMatrixTag(_assembly, _var.number(), _var.number());
134 :
135 458648 : precalculateJacobian();
136 2298424 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
137 : {
138 1839776 : initQpJacobian();
139 9201952 : for (_i = 0; _i < _test.size(); _i++)
140 37187200 : for (_j = 0; _j < _phi.size(); _j++)
141 : {
142 29825024 : _work_vector = computeQpJacobian() * _JxW[_qp] * _coord[_qp];
143 119300096 : _assembly.saveDiagLocalArrayJacobian(
144 29825024 : _local_ke, _i, _test.size(), _j, _phi.size(), _var.number(), _work_vector);
145 : }
146 : }
147 :
148 458648 : accumulateTaggedLocalMatrix();
149 :
150 458648 : if (_has_diag_save_in)
151 : {
152 0 : DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
153 0 : libMesh::Threads::spin_mutex::scoped_lock lock(libMesh::Threads::spin_mtx);
154 0 : for (const auto & var : _diag_save_in)
155 : {
156 0 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
157 0 : if (avar)
158 0 : avar->addSolution(diag);
159 : else
160 0 : mooseError("Save-in variable for an array kernel must be an array variable");
161 : }
162 0 : }
163 458648 : }
164 :
165 : RealEigenVector
166 21967648 : ArrayKernel::computeQpJacobian()
167 : {
168 43935296 : return RealEigenVector::Zero(_var.count());
169 : }
170 :
171 : void
172 650286 : ArrayKernel::computeOffDiagJacobian(const unsigned int jvar_num)
173 : {
174 650286 : const auto & jvar = getVariable(jvar_num);
175 :
176 650286 : bool same_var = (jvar_num == _var.number());
177 :
178 650286 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
179 :
180 : // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting on
181 : // the displaced mesh
182 650286 : auto phi_size = jvar.dofIndices().size();
183 :
184 650286 : precalculateOffDiagJacobian(jvar_num);
185 3290850 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
186 : {
187 2640564 : initQpOffDiagJacobian(jvar);
188 13544316 : for (_i = 0; _i < _test.size(); _i++)
189 57859560 : for (_j = 0; _j < phi_size; _j++)
190 : {
191 46955808 : _work_matrix = computeQpOffDiagJacobian(jvar) * _JxW[_qp] * _coord[_qp];
192 140867424 : _assembly.saveFullLocalArrayJacobian(
193 46955808 : _local_ke, _i, _test.size(), _j, phi_size, _var.number(), jvar_num, _work_matrix);
194 : }
195 : }
196 :
197 650286 : accumulateTaggedLocalMatrix();
198 :
199 650286 : if (_has_diag_save_in && same_var)
200 : {
201 256 : DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
202 256 : libMesh::Threads::spin_mutex::scoped_lock lock(libMesh::Threads::spin_mtx);
203 512 : for (const auto & var : _diag_save_in)
204 : {
205 256 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
206 256 : if (avar)
207 256 : avar->addSolution(diag);
208 : else
209 0 : mooseError("Save-in variable for an array kernel must be an array variable");
210 : }
211 256 : }
212 650286 : }
213 :
214 : RealEigenMatrix
215 44484392 : ArrayKernel::computeQpOffDiagJacobian(const MooseVariableFEBase & jvar)
216 : {
217 44484392 : if (jvar.number() == _var.number())
218 84967312 : return computeQpJacobian().asDiagonal();
219 : else
220 4001472 : return RealEigenMatrix::Zero(_var.count(), jvar.count());
221 : }
222 :
223 : void
224 0 : ArrayKernel::computeOffDiagJacobianScalar(unsigned int jvar)
225 : {
226 0 : MooseVariableScalar & jv = _sys.getScalarVariable(_tid, jvar);
227 0 : prepareMatrixTag(_assembly, _var.number(), jvar);
228 :
229 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
230 0 : for (_i = 0; _i < _test.size(); _i++)
231 : {
232 0 : _work_matrix = computeQpOffDiagJacobianScalar(jv) * _JxW[_qp] * _coord[_qp];
233 0 : _assembly.saveFullLocalArrayJacobian(
234 0 : _local_ke, _i, _test.size(), 0, 1, _var.number(), jvar, _work_matrix);
235 : }
236 :
237 0 : accumulateTaggedLocalMatrix();
238 0 : }
239 :
240 : RealEigenMatrix
241 0 : ArrayKernel::computeQpOffDiagJacobianScalar(const MooseVariableScalar & jvar)
242 : {
243 0 : return RealEigenMatrix::Zero(_var.count(), (unsigned int)jvar.order() + 1);
244 : }
|
