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 "DiffusionLHDGAssemblyHelper.h"
11 : #include "MooseVariableDependencyInterface.h"
12 : #include "MooseVariableFE.h"
13 : #include "SystemBase.h"
14 : #include "MooseObject.h"
15 : #include "MaterialPropertyInterface.h"
16 : #include "NonlinearThread.h"
17 : #include "TransientInterface.h"
18 :
19 : using namespace libMesh;
20 :
21 : InputParameters
22 9811 : DiffusionLHDGAssemblyHelper::validParams()
23 : {
24 9811 : auto params = emptyInputParameters();
25 39244 : params.addRequiredParam<NonlinearVariableName>(
26 : "gradient_variable", "The gradient of the diffusing specie concentration");
27 39244 : params.addRequiredParam<NonlinearVariableName>(
28 : "face_variable", "The concentration of the diffusing specie on faces");
29 39244 : params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
30 19622 : params.addParam<Real>("tau",
31 19622 : 1,
32 : "The stabilization coefficient required for discontinuous Galerkin "
33 : "schemes. This may be set to 0 for a mixed method with Raviart-Thomas.");
34 9811 : return params;
35 0 : }
36 :
37 316 : DiffusionLHDGAssemblyHelper::DiffusionLHDGAssemblyHelper(
38 : const MooseObject * const moose_obj,
39 : MaterialPropertyInterface * const mpi,
40 : MooseVariableDependencyInterface * const mvdi,
41 : const TransientInterface * const ti,
42 : const FEProblemBase & fe_problem,
43 : SystemBase & sys,
44 316 : const THREAD_ID tid)
45 : : ADFunctorInterface(moose_obj),
46 316 : _u_var(sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("variable"))),
47 632 : _grad_u_var(sys.getFieldVariable<RealVectorValue>(
48 632 : tid, moose_obj->getParam<NonlinearVariableName>("gradient_variable"))),
49 632 : _u_face_var(sys.getFieldVariable<Real>(
50 632 : tid, moose_obj->getParam<NonlinearVariableName>("face_variable"))),
51 316 : _qu_dof_indices(_grad_u_var.dofIndices()),
52 316 : _u_dof_indices(_u_var.dofIndices()),
53 316 : _lm_u_dof_indices(_u_face_var.dofIndices()),
54 316 : _qu_sol(_grad_u_var.sln()),
55 316 : _u_sol(_u_var.sln()),
56 316 : _lm_u_sol(_u_face_var.sln()),
57 316 : _vector_phi(_grad_u_var.phi()),
58 316 : _scalar_phi(_u_var.phi()),
59 316 : _grad_scalar_phi(_u_var.gradPhi()),
60 316 : _div_vector_phi(_grad_u_var.divPhi()),
61 316 : _vector_phi_face(_grad_u_var.phiFace()),
62 316 : _scalar_phi_face(_u_var.phiFace()),
63 316 : _lm_phi_face(_u_face_var.phiFace()),
64 632 : _diff(mpi->getMaterialProperty<Real>("diffusivity")),
65 316 : _ti(*ti),
66 632 : _tau(moose_obj->getParam<Real>("tau")),
67 316 : _cached_elem(nullptr),
68 316 : _moose_obj(*moose_obj),
69 316 : _dhah_fe_problem(fe_problem),
70 632 : _dhah_sys(sys)
71 : {
72 316 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<RealVectorValue> &>(_grad_u_var));
73 316 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_u_face_var));
74 316 : }
75 :
76 : void
77 316 : DiffusionLHDGAssemblyHelper::checkCoupling()
78 : {
79 : // This check must occur after FEProblemBase::init()
80 316 : if (_dhah_fe_problem.coupling() == Moose::COUPLING_FULL)
81 0 : return;
82 316 : else if (_dhah_fe_problem.coupling() == Moose::COUPLING_CUSTOM)
83 : {
84 316 : const auto * const cm = _dhah_fe_problem.couplingMatrix(_dhah_sys.number());
85 1264 : for (const auto i : make_range(cm->size()))
86 3792 : for (const auto j : make_range(cm->size()))
87 2844 : if ((*cm)(i, j) != true)
88 0 : goto error;
89 :
90 316 : return;
91 : }
92 :
93 0 : error:
94 0 : _moose_obj.mooseError(
95 : "This class encodes the full Jacobian regardless of user input file specification, "
96 : "so please request full coupling for system ",
97 0 : _dhah_sys.name(),
98 : " in your Preconditioning block for consistency");
99 : }
100 :
101 : void
102 156218 : DiffusionLHDGAssemblyHelper::vectorVolumeResidual(const MooseArray<Gradient> & vector_sol,
103 : const MooseArray<Number> & scalar_sol,
104 : const MooseArray<Real> & JxW,
105 : const QBase & qrule,
106 : DenseVector<Number> & vector_re)
107 : {
108 781090 : for (const auto qp : make_range(qrule.n_points()))
109 : {
110 624872 : const auto vector_qp_term = JxW[qp] * vector_sol[qp];
111 624872 : const auto scalar_qp_term = JxW[qp] * scalar_sol[qp];
112 4060184 : for (const auto i : index_range(vector_re))
113 : {
114 : // Vector equation dependence on vector dofs
115 3435312 : vector_re(i) += _vector_phi[i][qp] * vector_qp_term;
116 :
117 : // Vector equation dependence on scalar dofs
118 3435312 : vector_re(i) += _div_vector_phi[i][qp] * scalar_qp_term;
119 : }
120 : }
121 156218 : }
122 :
123 : void
124 102740 : DiffusionLHDGAssemblyHelper::vectorVolumeJacobian(const MooseArray<Real> & JxW,
125 : const QBase & qrule,
126 : DenseMatrix<Number> & vector_vector_jac,
127 : DenseMatrix<Number> & vector_scalar_jac)
128 : {
129 513700 : for (const auto qp : make_range(qrule.n_points()))
130 2661008 : for (const auto i : make_range(vector_vector_jac.m()))
131 : {
132 : // Vector equation dependence on vector dofs
133 2250048 : const auto vector_qpi_term = JxW[qp] * _vector_phi[i][qp];
134 16151136 : for (const auto j : make_range(vector_vector_jac.n()))
135 13901088 : vector_vector_jac(i, j) += vector_qpi_term * _vector_phi[j][qp];
136 :
137 : // Vector equation dependence on scalar dofs
138 2250048 : const auto scalar_qpi_term = JxW[qp] * _div_vector_phi[i][qp];
139 8568768 : for (const auto j : make_range(vector_scalar_jac.n()))
140 6318720 : vector_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi[j][qp];
141 : }
142 102740 : }
143 :
144 : void
145 156218 : DiffusionLHDGAssemblyHelper::scalarVolumeResidual(const MooseArray<Gradient> & vector_field,
146 : const Moose::Functor<Real> & source,
147 : const MooseArray<Real> & JxW,
148 : const QBase & qrule,
149 : const Elem * const current_elem,
150 : const MooseArray<Point> & q_point,
151 : DenseVector<Number> & scalar_re)
152 : {
153 781090 : for (const auto qp : make_range(qrule.n_points()))
154 : {
155 624872 : const auto vector_qp_term = JxW[qp] * _diff[qp] * vector_field[qp];
156 : // Evaluate source
157 : const auto f =
158 624872 : source(Moose::ElemQpArg{current_elem, qp, &qrule, q_point[qp]}, _ti.determineState());
159 624872 : const auto source_qp_term = JxW[qp] * f;
160 :
161 2093696 : for (const auto i : index_range(scalar_re))
162 : {
163 1468824 : scalar_re(i) += _grad_scalar_phi[i][qp] * vector_qp_term;
164 :
165 : // Scalar equation RHS
166 1468824 : scalar_re(i) -= _scalar_phi[i][qp] * source_qp_term;
167 : }
168 : }
169 156218 : }
170 :
171 : void
172 102740 : DiffusionLHDGAssemblyHelper::scalarVolumeJacobian(const MooseArray<Real> & JxW,
173 : const QBase & qrule,
174 : DenseMatrix<Number> & scalar_vector_jac)
175 : {
176 513700 : for (const auto qp : make_range(qrule.n_points()))
177 : {
178 410960 : const auto qp_term = JxW[qp] * _diff[qp];
179 1369040 : for (const auto i : make_range(scalar_vector_jac.m()))
180 : {
181 958080 : const auto qpi_term = qp_term * _grad_scalar_phi[i][qp];
182 : // Scalar equation dependence on vector dofs
183 7276800 : for (const auto j : make_range(scalar_vector_jac.n()))
184 6318720 : scalar_vector_jac(i, j) += qpi_term * _vector_phi[j][qp];
185 : }
186 : }
187 102740 : }
188 :
189 : void
190 551931 : DiffusionLHDGAssemblyHelper::vectorFaceResidual(const MooseArray<Number> & lm_sol,
191 : const MooseArray<Real> & JxW_face,
192 : const QBase & qrule_face,
193 : const MooseArray<Point> & normals,
194 : DenseVector<Number> & vector_re)
195 : {
196 : // Vector equation dependence on LM dofs
197 1655793 : for (const auto qp : make_range(qrule_face.n_points()))
198 : {
199 1103862 : const auto qp_term = JxW_face[qp] * lm_sol[qp] * normals[qp];
200 7270842 : for (const auto i : index_range(vector_re))
201 6166980 : vector_re(i) -= _vector_phi_face[i][qp] * qp_term;
202 : }
203 551931 : }
204 :
205 : void
206 364014 : DiffusionLHDGAssemblyHelper::vectorFaceJacobian(const MooseArray<Real> & JxW_face,
207 : const QBase & qrule_face,
208 : const MooseArray<Point> & normals,
209 : DenseMatrix<Number> & vector_lm_jac)
210 : {
211 1092042 : for (const auto qp : make_range(qrule_face.n_points()))
212 : {
213 728028 : const auto qp_term = JxW_face[qp] * normals[qp];
214 : // Vector equation dependence on LM dofs
215 4782852 : for (const auto i : make_range(vector_lm_jac.m()))
216 : {
217 4054824 : const auto qpi_term = qp_term * _vector_phi_face[i][qp];
218 34527384 : for (const auto j : make_range(vector_lm_jac.n()))
219 30472560 : vector_lm_jac(i, j) -= qpi_term * _lm_phi_face[j][qp];
220 : }
221 : }
222 364014 : }
223 :
224 : void
225 551931 : DiffusionLHDGAssemblyHelper::scalarFaceResidual(const MooseArray<Gradient> & vector_sol,
226 : const MooseArray<Number> & scalar_sol,
227 : const MooseArray<Number> & lm_sol,
228 : const MooseArray<Real> & JxW_face,
229 : const QBase & qrule_face,
230 : const MooseArray<Point> & normals,
231 : DenseVector<Number> & scalar_re)
232 : {
233 1655793 : for (const auto qp : make_range(qrule_face.n_points()))
234 : {
235 : // vector
236 1103862 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
237 : // stabilization term
238 1103862 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
239 : // scalar from stabilization term
240 1103862 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
241 : // lm from stabilization term
242 1103862 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
243 3727320 : for (const auto i : index_range(scalar_re))
244 2623458 : scalar_re(i) += _scalar_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
245 : }
246 551931 : }
247 :
248 : void
249 364014 : DiffusionLHDGAssemblyHelper::scalarFaceJacobian(const MooseArray<Real> & JxW_face,
250 : const QBase & qrule_face,
251 : const MooseArray<Point> & normals,
252 : DenseMatrix<Number> & scalar_vector_jac,
253 : DenseMatrix<Number> & scalar_scalar_jac,
254 : DenseMatrix<Number> & scalar_lm_jac)
255 : {
256 1092042 : for (const auto qp : make_range(qrule_face.n_points()))
257 : {
258 728028 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
259 728028 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
260 :
261 2447232 : for (const auto i : make_range(scalar_vector_jac.m()))
262 : {
263 1719204 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
264 13256988 : for (const auto j : make_range(scalar_vector_jac.n()))
265 11537784 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
266 :
267 1719204 : const auto scalar_qpi_term = stab_qp_term * _scalar_phi_face[i][qp];
268 7179888 : for (const auto j : make_range(scalar_scalar_jac.n()))
269 5460684 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
270 14594268 : for (const auto j : make_range(scalar_lm_jac.n()))
271 12875064 : scalar_lm_jac(i, j) -= scalar_qpi_term * _lm_phi_face[j][qp];
272 : }
273 : }
274 364014 : }
275 :
276 : void
277 551931 : DiffusionLHDGAssemblyHelper::lmFaceResidual(const MooseArray<Gradient> & vector_sol,
278 : const MooseArray<Number> & scalar_sol,
279 : const MooseArray<Number> & lm_sol,
280 : const MooseArray<Real> & JxW_face,
281 : const QBase & qrule_face,
282 : const MooseArray<Point> & normals,
283 : DenseVector<Number> & lm_re)
284 : {
285 1655793 : for (const auto qp : make_range(qrule_face.n_points()))
286 : {
287 : // vector
288 1103862 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
289 : // stabilization term
290 1103862 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
291 : // scalar from stabilization term
292 1103862 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
293 : // lm from stabilization term
294 1103862 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
295 9263034 : for (const auto i : index_range(lm_re))
296 8159172 : lm_re(i) += _lm_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
297 : }
298 551931 : }
299 :
300 : void
301 364014 : DiffusionLHDGAssemblyHelper::lmFaceJacobian(const MooseArray<Real> & JxW_face,
302 : const QBase & qrule_face,
303 : const MooseArray<Point> & normals,
304 : DenseMatrix<Number> & lm_vec_jac,
305 : DenseMatrix<Number> & lm_scalar_jac,
306 : DenseMatrix<Number> & lm_lm_jac)
307 : {
308 1092042 : for (const auto qp : make_range(qrule_face.n_points()))
309 : {
310 728028 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
311 728028 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
312 :
313 6120132 : for (const auto i : make_range(lm_vec_jac.m()))
314 : {
315 5392104 : const auto vector_qpi_term = vector_qp_term * _lm_phi_face[i][qp];
316 35864664 : for (const auto j : make_range(lm_vec_jac.n()))
317 30472560 : lm_vec_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
318 :
319 5392104 : const auto lm_qpi_term = stab_qp_term * _lm_phi_face[i][qp];
320 18267168 : for (const auto j : make_range(lm_scalar_jac.n()))
321 12875064 : lm_scalar_jac(i, j) += lm_qpi_term * _scalar_phi_face[j][qp];
322 45936216 : for (const auto j : make_range(lm_lm_jac.n()))
323 40544112 : lm_lm_jac(i, j) -= lm_qpi_term * _lm_phi_face[j][qp];
324 : }
325 : }
326 364014 : }
327 :
328 : void
329 14979 : DiffusionLHDGAssemblyHelper::vectorDirichletResidual(const Moose::Functor<Real> & dirichlet_value,
330 : const MooseArray<Real> & JxW_face,
331 : const QBase & qrule_face,
332 : const MooseArray<Point> & normals,
333 : const Elem * const current_elem,
334 : const unsigned int current_side,
335 : const MooseArray<Point> & q_point_face,
336 : DenseVector<Number> & vector_re)
337 : {
338 44937 : for (const auto qp : make_range(qrule_face.n_points()))
339 : {
340 89874 : const auto scalar_value = dirichlet_value(
341 29958 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
342 29958 : _ti.determineState());
343 29958 : const auto qp_term = JxW_face[qp] * normals[qp] * scalar_value;
344 :
345 : // External boundary -> Dirichlet faces -> Vector equation RHS
346 195018 : for (const auto i : index_range(_qu_dof_indices))
347 165060 : vector_re(i) -= qp_term * _vector_phi_face[i][qp];
348 : }
349 14979 : }
350 :
351 : void
352 14979 : DiffusionLHDGAssemblyHelper::scalarDirichletResidual(const MooseArray<Gradient> & vector_sol,
353 : const MooseArray<Number> & scalar_sol,
354 : const Moose::Functor<Real> & dirichlet_value,
355 : const MooseArray<Real> & JxW_face,
356 : const QBase & qrule_face,
357 : const MooseArray<Point> & normals,
358 : const Elem * const current_elem,
359 : const unsigned int current_side,
360 : const MooseArray<Point> & q_point_face,
361 : DenseVector<Number> & scalar_re)
362 : {
363 44937 : for (const auto qp : make_range(qrule_face.n_points()))
364 : {
365 89874 : const auto scalar_value = dirichlet_value(
366 29958 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
367 29958 : _ti.determineState());
368 29958 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
369 29958 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
370 29958 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
371 29958 : const auto lm_qp_term = stab_qp_term * scalar_value;
372 :
373 101016 : for (const auto i : index_range(_u_dof_indices))
374 71058 : scalar_re(i) += (scalar_qp_term - vector_qp_term - lm_qp_term) * _scalar_phi_face[i][qp];
375 : }
376 14979 : }
377 :
378 : void
379 9702 : DiffusionLHDGAssemblyHelper::scalarDirichletJacobian(const MooseArray<Real> & JxW_face,
380 : const QBase & qrule_face,
381 : const MooseArray<Point> & normals,
382 : DenseMatrix<Number> & scalar_vector_jac,
383 : DenseMatrix<Number> & scalar_scalar_jac)
384 : {
385 29106 : for (const auto qp : make_range(qrule_face.n_points()))
386 : {
387 19404 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
388 19404 : const auto scalar_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
389 64800 : for (const auto i : index_range(_u_dof_indices))
390 : {
391 45396 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
392 343548 : for (const auto j : index_range(_qu_dof_indices))
393 298152 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
394 :
395 45396 : const auto scalar_qpi_term = scalar_qp_term * _scalar_phi_face[i][qp];
396 186768 : for (const auto j : index_range(_u_dof_indices))
397 141372 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
398 : }
399 : }
400 9702 : }
401 :
402 : void
403 14979 : DiffusionLHDGAssemblyHelper::createIdentityResidual(const MooseArray<Real> & JxW,
404 : const QBase & qrule,
405 : const MooseArray<std::vector<Real>> & phi,
406 : const MooseArray<Number> & sol,
407 : DenseVector<Number> & re)
408 : {
409 44937 : for (const auto qp : make_range(qrule.n_points()))
410 : {
411 29958 : const auto qp_term = JxW[qp] * sol[qp];
412 245802 : for (const auto i : index_range(phi))
413 215844 : re(i) -= phi[i][qp] * qp_term;
414 : }
415 14979 : }
416 :
417 : void
418 9702 : DiffusionLHDGAssemblyHelper::createIdentityJacobian(const MooseArray<Real> & JxW,
419 : const QBase & qrule,
420 : const MooseArray<std::vector<Real>> & phi,
421 : DenseMatrix<Number> & ke)
422 : {
423 29106 : for (const auto qp : make_range(qrule.n_points()))
424 159828 : for (const auto i : index_range(phi))
425 : {
426 140424 : const auto qpi_term = JxW[qp] * phi[i][qp];
427 1174968 : for (const auto j : index_range(phi))
428 1034544 : ke(i, j) -= phi[j][qp] * qpi_term;
429 : }
430 9702 : }
|
