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 "ADKernel.h"
11 : #include "Assembly.h"
12 : #include "MooseVariable.h"
13 : #include "Problem.h"
14 : #include "SubProblem.h"
15 : #include "NonlinearSystemBase.h"
16 : #include "ADUtils.h"
17 : #include "FEProblemBase.h"
18 :
19 : // libmesh includes
20 : #include "libmesh/threads.h"
21 :
22 : template <typename T>
23 : InputParameters
24 493715 : ADKernelTempl<T>::validParams()
25 : {
26 493715 : auto params = KernelBase::validParams();
27 : if (std::is_same<T, Real>::value)
28 450521 : params.registerBase("Kernel");
29 : else if (std::is_same<T, RealVectorValue>::value)
30 43194 : params.registerBase("VectorKernel");
31 : else
32 : ::mooseError("unsupported ADKernelTempl specialization");
33 493715 : params += ADFunctorInterface::validParams();
34 493715 : return params;
35 0 : }
36 :
37 : template <typename T>
38 4535 : ADKernelTempl<T>::ADKernelTempl(const InputParameters & parameters)
39 : : KernelBase(parameters),
40 : MooseVariableInterface<T>(this,
41 : false,
42 : "variable",
43 : Moose::VarKindType::VAR_SOLVER,
44 : std::is_same<T, Real>::value ? Moose::VarFieldType::VAR_FIELD_STANDARD
45 : : Moose::VarFieldType::VAR_FIELD_VECTOR),
46 : ADFunctorInterface(this),
47 9070 : _var(*this->mooseVariable()),
48 4535 : _test(_var.phi()),
49 4535 : _grad_test(_var.adGradPhi()),
50 4535 : _regular_grad_test(_var.gradPhi()),
51 4535 : _u(_var.adSln()),
52 4535 : _grad_u(_var.adGradSln()),
53 4535 : _ad_JxW(_assembly.adJxW()),
54 4535 : _ad_coord(_assembly.adCoordTransformation()),
55 4535 : _ad_q_point(_assembly.adQPoints()),
56 4535 : _phi(_assembly.phi(_var)),
57 4535 : _grad_phi(_assembly.template adGradPhi<T>(_var)),
58 4535 : _regular_grad_phi(_assembly.gradPhi(_var)),
59 9070 : _my_elem(nullptr)
60 : {
61 4535 : _subproblem.haveADObjects(true);
62 :
63 4535 : addMooseVariableDependency(this->mooseVariable());
64 4535 : _save_in.resize(_save_in_strings.size());
65 4535 : _diag_save_in.resize(_diag_save_in_strings.size());
66 :
67 4535 : for (unsigned int i = 0; i < _save_in_strings.size(); i++)
68 : {
69 0 : MooseVariable * var = &_subproblem.getStandardVariable(_tid, _save_in_strings[i]);
70 :
71 0 : if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_save_in_strings[i]))
72 0 : paramError("save_in", "cannot use solution variable as save-in variable");
73 :
74 0 : if (var->feType() != _var.feType())
75 0 : paramError(
76 : "save_in",
77 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
78 0 : moose::internal::incompatVarMsg(*var, _var));
79 :
80 0 : _save_in[i] = var;
81 0 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
82 0 : addMooseVariableDependency(var);
83 : }
84 :
85 4535 : _has_save_in = _save_in.size() > 0;
86 :
87 4561 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
88 : {
89 26 : MooseVariable * var = &_subproblem.getStandardVariable(_tid, _diag_save_in_strings[i]);
90 :
91 26 : if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_diag_save_in_strings[i]))
92 0 : paramError("diag_save_in", "cannot use solution variable as diag save-in variable");
93 :
94 26 : if (var->feType() != _var.feType())
95 0 : paramError(
96 : "diag_save_in",
97 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
98 0 : moose::internal::incompatVarMsg(*var, _var));
99 :
100 26 : _diag_save_in[i] = var;
101 26 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
102 26 : addMooseVariableDependency(var);
103 : }
104 :
105 4535 : _has_diag_save_in = _diag_save_in.size() > 0;
106 4535 : }
107 :
108 : template <typename T>
109 : void
110 6374650 : ADKernelTempl<T>::computeResidual()
111 : {
112 6374650 : precalculateResidual();
113 :
114 6374650 : if (_residuals_nonad.size() != _test.size())
115 1590 : _residuals_nonad.resize(_test.size(), 0);
116 36359597 : for (auto & r : _residuals_nonad)
117 29984947 : r = 0;
118 :
119 34874769 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
120 : {
121 28500119 : const auto jxw_c = _JxW[_qp] * _coord[_qp];
122 180530970 : for (_i = 0; _i < _test.size(); _i++)
123 152030851 : _residuals_nonad[_i] += jxw_c * raw_value(computeQpResidual());
124 : }
125 :
126 6374650 : addResiduals(_assembly, _residuals_nonad, _var.dofIndices(), _var.scalingFactor());
127 :
128 6374650 : if (_has_save_in)
129 0 : for (unsigned int i = 0; i < _save_in.size(); i++)
130 0 : _save_in[i]->sys().solution().add_vector(_residuals_nonad.data(), _save_in[i]->dofIndices());
131 6374650 : }
132 :
133 : template <typename T>
134 : void
135 631513 : ADKernelTempl<T>::computeResidualsForJacobian()
136 : {
137 631513 : if (_residuals.size() != _test.size())
138 1713 : _residuals.resize(_test.size(), 0);
139 3921233 : for (auto & r : _residuals)
140 3289720 : r = 0;
141 :
142 631513 : precalculateResidual();
143 631513 : if (_use_displaced_mesh)
144 3500 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
145 : {
146 2800 : _r = _ad_JxW[_qp];
147 2800 : _r *= _ad_coord[_qp];
148 14000 : for (_i = 0; _i < _test.size(); _i++)
149 11200 : _residuals[_i] += _r * computeQpResidual();
150 : }
151 : else
152 3196325 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
153 : {
154 2565512 : const auto jxw_c = _JxW[_qp] * _coord[_qp];
155 16727252 : for (_i = 0; _i < _test.size(); _i++)
156 14161740 : _residuals[_i] += jxw_c * computeQpResidual();
157 : }
158 631513 : }
159 :
160 : template <typename T>
161 : void
162 2956 : ADKernelTempl<T>::computeJacobian()
163 : {
164 2956 : computeADJacobian();
165 :
166 2956 : if (_has_diag_save_in && !_sys.computingScalingJacobian())
167 0 : mooseError("_local_ke not computed for global AD indexing. Save-in is deprecated anyway. Use "
168 : "the tagging system instead.");
169 2956 : }
170 :
171 : template <typename T>
172 : void
173 1545319 : ADKernelTempl<T>::computeADJacobian()
174 : {
175 1545319 : computeResidualsForJacobian();
176 1545319 : addJacobian(_assembly, _residuals, dofIndices(), _var.scalingFactor());
177 1545319 : }
178 :
179 : template <typename T>
180 : void
181 25098 : ADKernelTempl<T>::jacobianSetup()
182 : {
183 25098 : _my_elem = nullptr;
184 25098 : }
185 :
186 : template <typename T>
187 : void
188 2218991 : ADKernelTempl<T>::computeOffDiagJacobian(const unsigned int)
189 : {
190 2218991 : if (_my_elem != _current_elem)
191 : {
192 1542363 : computeADJacobian();
193 1542363 : _my_elem = _current_elem;
194 : }
195 2218991 : }
196 :
197 : template <typename T>
198 : void
199 0 : ADKernelTempl<T>::computeOffDiagJacobianScalar(unsigned int /*jvar*/)
200 : {
201 0 : }
202 :
203 : template <typename T>
204 : void
205 3008 : ADKernelTempl<T>::computeResidualAndJacobian()
206 : {
207 3008 : computeResidualsForJacobian();
208 3008 : addResidualsAndJacobian(_assembly, _residuals, _var.dofIndices(), _var.scalingFactor());
209 3008 : }
210 :
211 : template class ADKernelTempl<Real>;
212 : template class ADKernelTempl<RealVectorValue>;
|
