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 "ADDGKernel.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 :
18 : // libmesh includes
19 : #include "libmesh/threads.h"
20 :
21 : InputParameters
22 9502 : ADDGKernel::validParams()
23 : {
24 9502 : InputParameters params = DGKernelBase::validParams();
25 9502 : params.addClassDescription(
26 : "Base class for all DG kernels making use of automatic differentiation");
27 9502 : return params;
28 0 : }
29 :
30 170 : ADDGKernel::ADDGKernel(const InputParameters & parameters)
31 : : DGKernelBase(parameters),
32 : NeighborMooseVariableInterface(
33 : this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
34 340 : _var(*mooseVariable()),
35 170 : _phi(_assembly.phiFace(_var)),
36 170 : _grad_phi(_assembly.gradPhiFace(_var)),
37 :
38 170 : _test(_var.phiFace()),
39 170 : _grad_test(_var.gradPhiFace()),
40 :
41 170 : _phi_neighbor(_assembly.phiFaceNeighbor(_var)),
42 170 : _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_var)),
43 :
44 170 : _test_neighbor(_var.phiFaceNeighbor()),
45 170 : _grad_test_neighbor(_var.gradPhiFaceNeighbor()),
46 :
47 170 : _u(_var.adSln()),
48 170 : _grad_u(_var.adGradSln()),
49 170 : _u_neighbor(_var.adSlnNeighbor()),
50 340 : _grad_u_neighbor(_var.adGradSlnNeighbor())
51 : {
52 170 : _subproblem.haveADObjects(true);
53 :
54 170 : addMooseVariableDependency(mooseVariable());
55 :
56 170 : _save_in.resize(_save_in_strings.size());
57 170 : _diag_save_in.resize(_diag_save_in_strings.size());
58 :
59 170 : for (unsigned int i = 0; i < _save_in_strings.size(); i++)
60 : {
61 0 : MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
62 0 : _save_in_strings[i],
63 : Moose::VarKindType::VAR_AUXILIARY,
64 : Moose::VarFieldType::VAR_FIELD_STANDARD);
65 :
66 0 : if (_sys.hasVariable(_save_in_strings[i]))
67 0 : mooseError("Trying to use solution variable " + _save_in_strings[i] +
68 0 : " as a save_in variable in " + name());
69 :
70 0 : if (var->feType() != _var.feType())
71 0 : paramError(
72 : "save_in",
73 : "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
74 0 : moose::internal::incompatVarMsg(*var, _var));
75 :
76 0 : _save_in[i] = var;
77 0 : var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
78 0 : addMooseVariableDependency(var);
79 : }
80 :
81 170 : _has_save_in = _save_in.size() > 0;
82 :
83 170 : for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
84 : {
85 0 : MooseVariableFEBase * var = &_subproblem.getVariable(_tid,
86 0 : _diag_save_in_strings[i],
87 : Moose::VarKindType::VAR_SOLVER,
88 : Moose::VarFieldType::VAR_FIELD_STANDARD);
89 :
90 0 : if (_sys.hasVariable(_diag_save_in_strings[i]))
91 0 : mooseError("Trying to use solution variable " + _diag_save_in_strings[i] +
92 0 : " as a diag_save_in variable in " + name());
93 :
94 0 : 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 0 : _diag_save_in[i] = var;
101 0 : var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
102 0 : addMooseVariableDependency(var);
103 : }
104 :
105 170 : _has_diag_save_in = _diag_save_in.size() > 0;
106 170 : }
107 :
108 : void
109 2005228 : ADDGKernel::computeElemNeighResidual(Moose::DGResidualType type)
110 : {
111 : bool is_elem;
112 2005228 : if (type == Moose::Element)
113 1002614 : is_elem = true;
114 : else
115 1002614 : is_elem = false;
116 :
117 2005228 : const VariableTestValue & test_space = is_elem ? _test : _test_neighbor;
118 :
119 2005228 : if (is_elem)
120 1002614 : prepareVectorTag(_assembly, _var.number());
121 : else
122 1002614 : prepareVectorTagNeighbor(_assembly, _var.number());
123 :
124 2005228 : precalculateResidual();
125 5957384 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
126 17917200 : for (_i = 0; _i < test_space.size(); _i++)
127 13965044 : _local_re(_i) += raw_value(_JxW[_qp] * _coord[_qp] * computeQpResidual(type));
128 :
129 2005228 : accumulateTaggedLocalResidual();
130 :
131 2005228 : if (_has_save_in)
132 0 : for (const auto & var : _save_in)
133 : {
134 : const std::vector<dof_id_type> & dof_indices =
135 0 : is_elem ? var->dofIndices() : var->dofIndicesNeighbor();
136 0 : var->sys().solution().add_vector(_local_re, dof_indices);
137 : }
138 2005228 : }
139 :
140 : void
141 0 : ADDGKernel::computeJacobian()
142 : {
143 : // AD only needs to do one computation for one variable because it does the derivatives all at
144 : // once
145 0 : if (!excludeBoundary())
146 : {
147 0 : computeElemNeighJacobian(Moose::ElementElement);
148 0 : computeElemNeighJacobian(Moose::NeighborNeighbor);
149 : }
150 0 : }
151 :
152 : void
153 0 : ADDGKernel::computeElemNeighJacobian(Moose::DGJacobianType type)
154 : {
155 : mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
156 : "With AD you should need one call per side");
157 :
158 0 : const VariableTestValue & test_space = type == Moose::ElementElement ? _test : _test_neighbor;
159 :
160 0 : std::vector<ADReal> residuals(test_space.size(), 0);
161 :
162 0 : precalculateResidual();
163 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
164 0 : for (_i = 0; _i < test_space.size(); _i++)
165 0 : residuals[_i] +=
166 0 : _JxW[_qp] * _coord[_qp] *
167 0 : computeQpResidual(type == Moose::ElementElement ? Moose::Element : Moose::Neighbor);
168 :
169 0 : addJacobian(_assembly,
170 : residuals,
171 0 : type == Moose::ElementElement ? _var.dofIndices() : _var.dofIndicesNeighbor(),
172 0 : _var.scalingFactor());
173 :
174 0 : if (_has_diag_save_in)
175 : {
176 0 : unsigned int rows = _local_ke.m();
177 0 : DenseVector<Number> diag(rows);
178 0 : for (unsigned int i = 0; i < rows; i++)
179 0 : diag(i) = _local_ke(i, i);
180 :
181 0 : for (const auto & var : _diag_save_in)
182 : {
183 0 : if (type == Moose::ElementElement)
184 0 : var->sys().solution().add_vector(diag, var->dofIndices());
185 : else
186 0 : var->sys().solution().add_vector(diag, var->dofIndicesNeighbor());
187 : }
188 0 : }
189 0 : }
190 :
191 : void
192 40840 : ADDGKernel::computeOffDiagJacobian(const unsigned int jvar_num)
193 : {
194 : // AD only needs to do one computation for one variable because it does the derivatives all at
195 : // once
196 40840 : if (!excludeBoundary() && jvar_num == _var.number())
197 : {
198 34362 : const auto & jvar = getVariable(jvar_num);
199 34362 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
200 34362 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
201 : }
202 40840 : }
203 :
204 : void
205 68724 : ADDGKernel::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, const MooseVariableFEBase &)
206 : {
207 : mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
208 : "With AD you should need one call per side");
209 :
210 68724 : const VariableTestValue & test_space = type == Moose::ElementElement ? _test : _test_neighbor;
211 :
212 68724 : std::vector<ADReal> residuals(test_space.size(), 0);
213 :
214 68724 : precalculateResidual();
215 203176 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
216 598088 : for (_i = 0; _i < test_space.size(); _i++)
217 463636 : residuals[_i] +=
218 927272 : _JxW[_qp] * _coord[_qp] *
219 1390908 : computeQpResidual(type == Moose::ElementElement ? Moose::Element : Moose::Neighbor);
220 :
221 137448 : addJacobian(_assembly,
222 : residuals,
223 68724 : type == Moose::ElementElement ? _var.dofIndices() : _var.dofIndicesNeighbor(),
224 68724 : _var.scalingFactor());
225 68724 : }
|
