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 494536 : ADKernelTempl<T>::validParams()
25 : {
26 494536 : auto params = KernelBase::validParams();
27 : if (std::is_same<T, Real>::value)
28 902644 : params.registerBase("Kernel");
29 : else if (std::is_same<T, RealVectorValue>::value)
30 86428 : params.registerBase("VectorKernel");
31 : else
32 : ::mooseError("unsupported ADKernelTempl specialization");
33 494536 : params += ADFunctorInterface::validParams();
34 494536 : return params;
35 0 : }
36 :
37 : template <typename T>
38 4955 : 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 9910 : _var(*this->mooseVariable()),
48 4955 : _test(_var.phi()),
49 4955 : _grad_test(_var.adGradPhi()),
50 4955 : _regular_grad_test(_var.gradPhi()),
51 4955 : _u(_var.adSln()),
52 4955 : _grad_u(_var.adGradSln()),
53 4955 : _ad_JxW(_assembly.adJxW()),
54 4955 : _ad_coord(_assembly.adCoordTransformation()),
55 4955 : _ad_q_point(_assembly.adQPoints()),
56 4955 : _phi(_assembly.phi(_var)),
57 4955 : _grad_phi(_assembly.template adGradPhi<T>(_var)),
58 4955 : _regular_grad_phi(_assembly.gradPhi(_var)),
59 19820 : _my_elem(nullptr)
60 : {
61 4955 : _subproblem.haveADObjects(true);
62 :
63 4955 : addMooseVariableDependency(this->mooseVariable());
64 4955 : _save_in.resize(_save_in_strings.size());
65 4955 : _diag_save_in.resize(_diag_save_in_strings.size());
66 :
67 4955 : 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 4955 : _has_save_in = _save_in.size() > 0;
86 :
87 4983 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
88 : {
89 28 : MooseVariable * var = &_subproblem.getStandardVariable(_tid, _diag_save_in_strings[i]);
90 :
91 28 : 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 28 : 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 28 : _diag_save_in[i] = var;
101 28 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
102 28 : addMooseVariableDependency(var);
103 : }
104 :
105 4955 : _has_diag_save_in = _diag_save_in.size() > 0;
106 4955 : }
107 :
108 : template <typename T>
109 : void
110 6599714 : ADKernelTempl<T>::computeResidual()
111 : {
112 6599714 : precalculateResidual();
113 :
114 6599714 : if (_residuals_nonad.size() != _test.size())
115 1676 : _residuals_nonad.resize(_test.size(), 0);
116 37698657 : for (auto & r : _residuals_nonad)
117 31098943 : r = 0;
118 :
119 36065394 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
120 : {
121 29465680 : const auto jxw_c = _JxW[_qp] * _coord[_qp];
122 186598675 : for (_i = 0; _i < _test.size(); _i++)
123 157132995 : _residuals_nonad[_i] += jxw_c * raw_value(computeQpResidual());
124 : }
125 :
126 6599714 : addResiduals(_assembly, _residuals_nonad, _var.dofIndices(), _var.scalingFactor());
127 :
128 6599714 : 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 6599714 : }
132 :
133 : template <typename T>
134 : void
135 713413 : ADKernelTempl<T>::computeResidualsForJacobian()
136 : {
137 713413 : if (_residuals.size() != _test.size())
138 1802 : _residuals.resize(_test.size(), 0);
139 4437453 : for (auto & r : _residuals)
140 3724040 : r = 0;
141 :
142 713413 : precalculateResidual();
143 713413 : 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 3611035 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
153 : {
154 2898322 : const auto jxw_c = _JxW[_qp] * _coord[_qp];
155 18924982 : for (_i = 0; _i < _test.size(); _i++)
156 16026660 : _residuals[_i] += jxw_c * computeQpResidual();
157 : }
158 713413 : }
159 :
160 : template <typename T>
161 : void
162 3164 : ADKernelTempl<T>::computeJacobian()
163 : {
164 3164 : computeADJacobian();
165 :
166 3164 : 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 3164 : }
170 :
171 : template <typename T>
172 : void
173 1759100 : ADKernelTempl<T>::computeADJacobian()
174 : {
175 1759100 : computeResidualsForJacobian();
176 1759100 : addJacobian(_assembly, _residuals, dofIndices(), _var.scalingFactor());
177 1759100 : }
178 :
179 : template <typename T>
180 : void
181 28495 : ADKernelTempl<T>::jacobianSetup()
182 : {
183 28495 : _my_elem = nullptr;
184 28495 : }
185 :
186 : template <typename T>
187 : void
188 2520191 : ADKernelTempl<T>::computeOffDiagJacobian(const unsigned int)
189 : {
190 2520191 : if (_my_elem != _current_elem)
191 : {
192 1755936 : computeADJacobian();
193 1755936 : _my_elem = _current_elem;
194 : }
195 2520191 : }
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 3048 : ADKernelTempl<T>::computeResidualAndJacobian()
206 : {
207 3048 : computeResidualsForJacobian();
208 3048 : addResidualsAndJacobian(_assembly, _residuals, _var.dofIndices(), _var.scalingFactor());
209 3048 : }
210 :
211 : template class ADKernelTempl<Real>;
212 : template class ADKernelTempl<RealVectorValue>;
|
