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 103666 : ArrayKernel::validParams()
24 : {
25 103666 : InputParameters params = KernelBase::validParams();
26 103666 : params.registerBase("ArrayKernel");
27 103666 : return params;
28 0 : }
29 :
30 1953 : 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 3906 : _var(*mooseVariable()),
38 1953 : _test(_var.phi()),
39 1953 : _grad_test(_var.gradPhi()),
40 1953 : _array_grad_test(_var.arrayGradPhi()),
41 1953 : _phi(_assembly.phi(_var)),
42 1953 : _grad_phi(_assembly.gradPhi(_var)),
43 1953 : _u(_is_implicit ? _var.sln() : _var.slnOld()),
44 1953 : _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
45 1953 : _count(_var.count()),
46 3906 : _work_vector(_count)
47 : {
48 1953 : addMooseVariableDependency(mooseVariable());
49 :
50 1953 : _save_in.resize(_save_in_strings.size());
51 1953 : _diag_save_in.resize(_diag_save_in_strings.size());
52 :
53 1967 : for (unsigned int i = 0; i < _save_in_strings.size(); i++)
54 : {
55 14 : ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _save_in_strings[i]);
56 :
57 14 : 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 14 : 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 14 : _save_in[i] = var;
67 14 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
68 14 : addMooseVariableDependency(var);
69 : }
70 :
71 1953 : _has_save_in = _save_in.size() > 0;
72 :
73 1967 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
74 : {
75 14 : ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _diag_save_in_strings[i]);
76 :
77 14 : 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 14 : 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 14 : _diag_save_in[i] = var;
87 14 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
88 14 : addMooseVariableDependency(var);
89 : }
90 :
91 1953 : _has_diag_save_in = _diag_save_in.size() > 0;
92 1953 : }
93 :
94 : void
95 5157197 : ArrayKernel::computeResidual()
96 : {
97 5157197 : prepareVectorTag(_assembly, _var.number());
98 :
99 5157197 : precalculateResidual();
100 28869469 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
101 : {
102 23712272 : initQpResidual();
103 142947652 : for (_i = 0; _i < _test.size(); _i++)
104 : {
105 119235380 : _work_vector.setZero();
106 119235380 : 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 119235380 : _work_vector *= _JxW[_qp] * _coord[_qp];
110 119235380 : _assembly.saveLocalArrayResidual(_local_re, _i, _test.size(), _work_vector);
111 : }
112 : }
113 :
114 5157197 : accumulateTaggedLocalResidual();
115 :
116 5157197 : if (_has_save_in)
117 : {
118 432 : libMesh::Threads::spin_mutex::scoped_lock lock(libMesh::Threads::spin_mtx);
119 864 : for (const auto & var : _save_in)
120 : {
121 432 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
122 432 : if (avar)
123 432 : avar->addSolution(_local_re);
124 : else
125 0 : mooseError("Save-in variable for an array kernel must be an array variable");
126 : }
127 432 : }
128 5157197 : }
129 :
130 : void
131 529411 : ArrayKernel::computeJacobian()
132 : {
133 529411 : prepareMatrixTag(_assembly, _var.number(), _var.number());
134 :
135 529411 : precalculateJacobian();
136 2652887 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
137 : {
138 2123476 : initQpJacobian();
139 10618436 : for (_i = 0; _i < _test.size(); _i++)
140 42895760 : for (_j = 0; _j < _phi.size(); _j++)
141 : {
142 34400800 : _work_vector = computeQpJacobian() * _JxW[_qp] * _coord[_qp];
143 137603200 : _assembly.saveDiagLocalArrayJacobian(
144 34400800 : _local_ke, _i, _test.size(), _j, _phi.size(), _var.number(), _work_vector);
145 : }
146 : }
147 :
148 529411 : accumulateTaggedLocalMatrix();
149 :
150 529411 : 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 529411 : }
164 :
165 : RealEigenVector
166 25315732 : ArrayKernel::computeQpJacobian()
167 : {
168 50631464 : return RealEigenVector::Zero(_var.count());
169 : }
170 :
171 : void
172 744798 : ArrayKernel::computeOffDiagJacobian(const unsigned int jvar_num)
173 : {
174 744798 : const auto & jvar = getVariable(jvar_num);
175 :
176 744798 : bool same_var = (jvar_num == _var.number());
177 :
178 744798 : 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 744798 : auto phi_size = jvar.dofIndices().size();
183 :
184 744798 : precalculateOffDiagJacobian(jvar_num);
185 3765600 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
186 : {
187 3020802 : initQpOffDiagJacobian(jvar);
188 15465708 : for (_i = 0; _i < _test.size(); _i++)
189 65801880 : for (_j = 0; _j < phi_size; _j++)
190 : {
191 53356974 : _work_matrix = computeQpOffDiagJacobian(jvar) * _JxW[_qp] * _coord[_qp];
192 160070922 : _assembly.saveFullLocalArrayJacobian(
193 53356974 : _local_ke, _i, _test.size(), _j, phi_size, _var.number(), jvar_num, _work_matrix);
194 : }
195 : }
196 :
197 744798 : accumulateTaggedLocalMatrix();
198 :
199 744798 : if (_has_diag_save_in && same_var)
200 : {
201 288 : DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
202 288 : libMesh::Threads::spin_mutex::scoped_lock lock(libMesh::Threads::spin_mtx);
203 576 : for (const auto & var : _diag_save_in)
204 : {
205 288 : auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
206 288 : if (avar)
207 288 : avar->addSolution(diag);
208 : else
209 0 : mooseError("Save-in variable for an array kernel must be an array variable");
210 : }
211 288 : }
212 744798 : }
213 :
214 : RealEigenMatrix
215 50653845 : ArrayKernel::computeQpOffDiagJacobian(const MooseVariableFEBase & jvar)
216 : {
217 50653845 : if (jvar.number() == _var.number())
218 97166090 : return computeQpJacobian().asDiagonal();
219 : else
220 4141600 : 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 : }
|
