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 43042 : ADDGKernel::validParams()
23 : {
24 43042 : InputParameters params = DGKernelBase::validParams();
25 43042 : params.addClassDescription(
26 : "Base class for all DG kernels making use of automatic differentiation");
27 43042 : return params;
28 0 : }
29 :
30 131 : ADDGKernel::ADDGKernel(const InputParameters & parameters)
31 : : DGKernelBase(parameters),
32 : NeighborMooseVariableInterface(
33 : this, false, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
34 262 : _var(*mooseVariable()),
35 131 : _phi(_assembly.phiFace(_var)),
36 131 : _grad_phi(_assembly.gradPhiFace(_var)),
37 :
38 131 : _test(_var.phiFace()),
39 131 : _grad_test(_var.gradPhiFace()),
40 :
41 131 : _phi_neighbor(_assembly.phiFaceNeighbor(_var)),
42 131 : _grad_phi_neighbor(_assembly.gradPhiFaceNeighbor(_var)),
43 :
44 131 : _test_neighbor(_var.phiFaceNeighbor()),
45 131 : _grad_test_neighbor(_var.gradPhiFaceNeighbor()),
46 :
47 131 : _u(_var.adSln()),
48 131 : _grad_u(_var.adGradSln()),
49 131 : _u_neighbor(_var.adSlnNeighbor()),
50 262 : _grad_u_neighbor(_var.adGradSlnNeighbor())
51 : {
52 131 : _subproblem.haveADObjects(true);
53 :
54 131 : addMooseVariableDependency(mooseVariable());
55 :
56 131 : _save_in.resize(_save_in_strings.size());
57 131 : _diag_save_in.resize(_diag_save_in_strings.size());
58 :
59 131 : 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 131 : _has_save_in = _save_in.size() > 0;
82 :
83 131 : 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 131 : _has_diag_save_in = _diag_save_in.size() > 0;
106 131 : }
107 :
108 : void
109 1377736 : ADDGKernel::computeElemNeighResidual(Moose::DGResidualType type)
110 : {
111 : bool is_elem;
112 1377736 : if (type == Moose::Element)
113 688868 : is_elem = true;
114 : else
115 688868 : is_elem = false;
116 :
117 1377736 : const VariableTestValue & test_space = is_elem ? _test : _test_neighbor;
118 :
119 1377736 : if (is_elem)
120 688868 : prepareVectorTag(_assembly, _var.number());
121 : else
122 688868 : prepareVectorTagNeighbor(_assembly, _var.number());
123 :
124 4065696 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
125 13196324 : for (_i = 0; _i < test_space.size(); _i++)
126 10508364 : _local_re(_i) += raw_value(_JxW[_qp] * _coord[_qp] * computeQpResidual(type));
127 :
128 1377736 : accumulateTaggedLocalResidual();
129 :
130 1377736 : if (_has_save_in)
131 : {
132 0 : Threads::spin_mutex::scoped_lock lock(_resid_vars_mutex);
133 0 : for (const auto & var : _save_in)
134 : {
135 : const std::vector<dof_id_type> & dof_indices =
136 0 : is_elem ? var->dofIndices() : var->dofIndicesNeighbor();
137 0 : var->sys().solution().add_vector(_local_re, dof_indices);
138 : }
139 0 : }
140 1377736 : }
141 :
142 : void
143 0 : ADDGKernel::computeJacobian()
144 : {
145 : // AD only needs to do one computation for one variable because it does the derivatives all at
146 : // once
147 0 : if (!excludeBoundary())
148 : {
149 0 : computeElemNeighJacobian(Moose::ElementElement);
150 0 : computeElemNeighJacobian(Moose::NeighborNeighbor);
151 : }
152 0 : }
153 :
154 : void
155 0 : ADDGKernel::computeElemNeighJacobian(Moose::DGJacobianType type)
156 : {
157 : mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
158 : "With AD you should need one call per side");
159 :
160 0 : const VariableTestValue & test_space = type == Moose::ElementElement ? _test : _test_neighbor;
161 :
162 0 : std::vector<ADReal> residuals(test_space.size(), 0);
163 :
164 0 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
165 0 : for (_i = 0; _i < test_space.size(); _i++)
166 0 : residuals[_i] +=
167 0 : _JxW[_qp] * _coord[_qp] *
168 0 : computeQpResidual(type == Moose::ElementElement ? Moose::Element : Moose::Neighbor);
169 :
170 0 : addJacobian(_assembly,
171 : residuals,
172 0 : type == Moose::ElementElement ? _var.dofIndices() : _var.dofIndicesNeighbor(),
173 0 : _var.scalingFactor());
174 :
175 0 : if (_has_diag_save_in)
176 : {
177 0 : unsigned int rows = _local_ke.m();
178 0 : DenseVector<Number> diag(rows);
179 0 : for (unsigned int i = 0; i < rows; i++)
180 0 : diag(i) = _local_ke(i, i);
181 :
182 0 : Threads::spin_mutex::scoped_lock lock(_jacoby_vars_mutex);
183 0 : for (const auto & var : _diag_save_in)
184 : {
185 0 : if (type == Moose::ElementElement)
186 0 : var->sys().solution().add_vector(diag, var->dofIndices());
187 : else
188 0 : var->sys().solution().add_vector(diag, var->dofIndicesNeighbor());
189 : }
190 0 : }
191 0 : }
192 :
193 : void
194 32518 : ADDGKernel::computeOffDiagJacobian(const unsigned int jvar_num)
195 : {
196 : // AD only needs to do one computation for one variable because it does the derivatives all at
197 : // once
198 32518 : if (!excludeBoundary() && jvar_num == _var.number())
199 : {
200 26396 : const auto & jvar = getVariable(jvar_num);
201 26396 : computeOffDiagElemNeighJacobian(Moose::ElementElement, jvar);
202 26396 : computeOffDiagElemNeighJacobian(Moose::NeighborNeighbor, jvar);
203 : }
204 32518 : }
205 :
206 : void
207 52792 : ADDGKernel::computeOffDiagElemNeighJacobian(Moose::DGJacobianType type, const MooseVariableFEBase &)
208 : {
209 : mooseAssert(type == Moose::ElementElement || type == Moose::NeighborNeighbor,
210 : "With AD you should need one call per side");
211 :
212 52792 : const VariableTestValue & test_space = type == Moose::ElementElement ? _test : _test_neighbor;
213 :
214 52792 : std::vector<ADReal> residuals(test_space.size(), 0);
215 :
216 154992 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
217 468312 : for (_i = 0; _i < test_space.size(); _i++)
218 366112 : residuals[_i] +=
219 732224 : _JxW[_qp] * _coord[_qp] *
220 1098336 : computeQpResidual(type == Moose::ElementElement ? Moose::Element : Moose::Neighbor);
221 :
222 105584 : addJacobian(_assembly,
223 : residuals,
224 52792 : type == Moose::ElementElement ? _var.dofIndices() : _var.dofIndicesNeighbor(),
225 52792 : _var.scalingFactor());
226 52792 : }
|
