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 43831 : DiffusionLHDGAssemblyHelper::validParams()
23 : {
24 43831 : auto params = emptyInputParameters();
25 43831 : params.addRequiredParam<NonlinearVariableName>(
26 : "gradient_variable", "The gradient of the diffusing specie concentration");
27 43831 : params.addRequiredParam<NonlinearVariableName>(
28 : "face_variable", "The concentration of the diffusing specie on faces");
29 43831 : params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
30 131493 : params.addParam<Real>("tau",
31 87662 : 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 43831 : return params;
35 0 : }
36 :
37 520 : 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 520 : const THREAD_ID tid)
45 : : ADFunctorInterface(moose_obj),
46 520 : _u_var(sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("variable"))),
47 1040 : _grad_u_var(sys.getFieldVariable<RealVectorValue>(
48 520 : tid, moose_obj->getParam<NonlinearVariableName>("gradient_variable"))),
49 1040 : _u_face_var(sys.getFieldVariable<Real>(
50 520 : tid, moose_obj->getParam<NonlinearVariableName>("face_variable"))),
51 520 : _qu_dof_indices(_grad_u_var.dofIndices()),
52 520 : _u_dof_indices(_u_var.dofIndices()),
53 520 : _lm_u_dof_indices(_u_face_var.dofIndices()),
54 520 : _qu_sol(_grad_u_var.sln()),
55 520 : _u_sol(_u_var.sln()),
56 520 : _lm_u_sol(_u_face_var.sln()),
57 520 : _vector_phi(_grad_u_var.phi()),
58 520 : _scalar_phi(_u_var.phi()),
59 520 : _grad_scalar_phi(_u_var.gradPhi()),
60 520 : _div_vector_phi(_grad_u_var.divPhi()),
61 520 : _vector_phi_face(_grad_u_var.phiFace()),
62 520 : _scalar_phi_face(_u_var.phiFace()),
63 520 : _lm_phi_face(_u_face_var.phiFace()),
64 520 : _diff(mpi->getMaterialProperty<Real>("diffusivity")),
65 520 : _ti(*ti),
66 520 : _tau(moose_obj->getParam<Real>("tau")),
67 520 : _cached_elem(nullptr),
68 520 : _moose_obj(*moose_obj),
69 520 : _dhah_fe_problem(fe_problem),
70 1040 : _dhah_sys(sys)
71 : {
72 520 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<RealVectorValue> &>(_grad_u_var));
73 520 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_u_face_var));
74 520 : }
75 :
76 : void
77 520 : DiffusionLHDGAssemblyHelper::checkCoupling()
78 : {
79 : // This check must occur after FEProblemBase::init()
80 520 : if (_dhah_fe_problem.coupling() == Moose::COUPLING_FULL)
81 0 : return;
82 520 : else if (_dhah_fe_problem.coupling() == Moose::COUPLING_CUSTOM)
83 : {
84 520 : const auto * const cm = _dhah_fe_problem.couplingMatrix(_dhah_sys.number());
85 2080 : for (const auto i : make_range(cm->size()))
86 6240 : for (const auto j : make_range(cm->size()))
87 4680 : if ((*cm)(i, j) != true)
88 0 : goto error;
89 :
90 520 : 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 356368 : 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 1781840 : for (const auto qp : make_range(qrule.n_points()))
109 : {
110 1425472 : const auto vector_qp_term = JxW[qp] * vector_sol[qp];
111 1425472 : const auto scalar_qp_term = JxW[qp] * scalar_sol[qp];
112 9248512 : for (const auto i : index_range(vector_re))
113 : {
114 : // Vector equation dependence on vector dofs
115 7823040 : vector_re(i) += _vector_phi[i][qp] * vector_qp_term;
116 :
117 : // Vector equation dependence on scalar dofs
118 7823040 : vector_re(i) += _div_vector_phi[i][qp] * scalar_qp_term;
119 : }
120 : }
121 356368 : }
122 :
123 : void
124 236336 : DiffusionLHDGAssemblyHelper::vectorVolumeJacobian(const MooseArray<Real> & JxW,
125 : const QBase & qrule,
126 : DenseMatrix<Number> & vector_vector_jac,
127 : DenseMatrix<Number> & vector_scalar_jac)
128 : {
129 1181680 : for (const auto qp : make_range(qrule.n_points()))
130 6116768 : for (const auto i : make_range(vector_vector_jac.m()))
131 : {
132 : // Vector equation dependence on vector dofs
133 5171424 : const auto vector_qpi_term = JxW[qp] * _vector_phi[i][qp];
134 37143232 : for (const auto j : make_range(vector_vector_jac.n()))
135 31971808 : vector_vector_jac(i, j) += vector_qpi_term * _vector_phi[j][qp];
136 :
137 : // Vector equation dependence on scalar dofs
138 5171424 : const auto scalar_qpi_term = JxW[qp] * _div_vector_phi[i][qp];
139 19684416 : for (const auto j : make_range(vector_scalar_jac.n()))
140 14512992 : vector_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi[j][qp];
141 : }
142 236336 : }
143 :
144 : void
145 356368 : 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 1781840 : for (const auto qp : make_range(qrule.n_points()))
154 : {
155 1425472 : const auto vector_qp_term = JxW[qp] * _diff[qp] * vector_field[qp];
156 : // Evaluate source
157 : const auto f =
158 1425472 : source(Moose::ElemQpArg{current_elem, qp, &qrule, q_point[qp]}, _ti.determineState());
159 1425472 : const auto source_qp_term = JxW[qp] * f;
160 :
161 4756832 : for (const auto i : index_range(scalar_re))
162 : {
163 3331360 : scalar_re(i) += _grad_scalar_phi[i][qp] * vector_qp_term;
164 :
165 : // Scalar equation RHS
166 3331360 : scalar_re(i) -= _scalar_phi[i][qp] * source_qp_term;
167 : }
168 : }
169 356368 : }
170 :
171 : void
172 236336 : DiffusionLHDGAssemblyHelper::scalarVolumeJacobian(const MooseArray<Real> & JxW,
173 : const QBase & qrule,
174 : DenseMatrix<Number> & scalar_vector_jac)
175 : {
176 1181680 : for (const auto qp : make_range(qrule.n_points()))
177 : {
178 945344 : const auto qp_term = JxW[qp] * _diff[qp];
179 3141968 : for (const auto i : make_range(scalar_vector_jac.m()))
180 : {
181 2196624 : const auto qpi_term = qp_term * _grad_scalar_phi[i][qp];
182 : // Scalar equation dependence on vector dofs
183 16709616 : for (const auto j : make_range(scalar_vector_jac.n()))
184 14512992 : scalar_vector_jac(i, j) += qpi_term * _vector_phi[j][qp];
185 : }
186 : }
187 236336 : }
188 :
189 : void
190 1264929 : 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 3794787 : for (const auto qp : make_range(qrule_face.n_points()))
198 : {
199 2529858 : const auto qp_term = JxW_face[qp] * lm_sol[qp] * normals[qp];
200 16647918 : for (const auto i : index_range(vector_re))
201 14118060 : vector_re(i) -= _vector_phi_face[i][qp] * qp_term;
202 : }
203 1264929 : }
204 :
205 : void
206 839866 : 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 2519598 : for (const auto qp : make_range(qrule_face.n_points()))
212 : {
213 1679732 : const auto qp_term = JxW_face[qp] * normals[qp];
214 : // Vector equation dependence on LM dofs
215 11030348 : for (const auto i : make_range(vector_lm_jac.m()))
216 : {
217 9350616 : const auto qpi_term = qp_term * _vector_phi_face[i][qp];
218 79789416 : for (const auto j : make_range(vector_lm_jac.n()))
219 70438800 : vector_lm_jac(i, j) -= qpi_term * _lm_phi_face[j][qp];
220 : }
221 : }
222 839866 : }
223 :
224 : void
225 1264929 : 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 3794787 : for (const auto qp : make_range(qrule_face.n_points()))
234 : {
235 : // vector
236 2529858 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
237 : // stabilization term
238 2529858 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
239 : // scalar from stabilization term
240 2529858 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
241 : // lm from stabilization term
242 2529858 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
243 8516040 : for (const auto i : index_range(scalar_re))
244 5986182 : scalar_re(i) += _scalar_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
245 : }
246 1264929 : }
247 :
248 : void
249 839866 : 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 2519598 : for (const auto qp : make_range(qrule_face.n_points()))
257 : {
258 1679732 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
259 1679732 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
260 :
261 5636448 : for (const auto i : make_range(scalar_vector_jac.m()))
262 : {
263 3956716 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
264 30553012 : for (const auto j : make_range(scalar_vector_jac.n()))
265 26596296 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
266 :
267 3956716 : const auto scalar_qpi_term = stab_qp_term * _scalar_phi_face[i][qp];
268 16536272 : for (const auto j : make_range(scalar_scalar_jac.n()))
269 12579556 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
270 33669972 : for (const auto j : make_range(scalar_lm_jac.n()))
271 29713256 : scalar_lm_jac(i, j) -= scalar_qpi_term * _lm_phi_face[j][qp];
272 : }
273 : }
274 839866 : }
275 :
276 : void
277 1264929 : 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 3794787 : for (const auto qp : make_range(qrule_face.n_points()))
286 : {
287 : // vector
288 2529858 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
289 : // stabilization term
290 2529858 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
291 : // scalar from stabilization term
292 2529858 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
293 : // lm from stabilization term
294 2529858 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
295 21293006 : for (const auto i : index_range(lm_re))
296 18763148 : lm_re(i) += _lm_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
297 : }
298 1264929 : }
299 :
300 : void
301 839866 : 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 2519598 : for (const auto qp : make_range(qrule_face.n_points()))
309 : {
310 1679732 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
311 1679732 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
312 :
313 14147308 : for (const auto i : make_range(lm_vec_jac.m()))
314 : {
315 12467576 : const auto vector_qpi_term = vector_qp_term * _lm_phi_face[i][qp];
316 82906376 : for (const auto j : make_range(lm_vec_jac.n()))
317 70438800 : lm_vec_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
318 :
319 12467576 : const auto lm_qpi_term = stab_qp_term * _lm_phi_face[i][qp];
320 42180832 : for (const auto j : make_range(lm_scalar_jac.n()))
321 29713256 : lm_scalar_jac(i, j) += lm_qpi_term * _scalar_phi_face[j][qp];
322 106386504 : for (const auto j : make_range(lm_lm_jac.n()))
323 93918928 : lm_lm_jac(i, j) -= lm_qpi_term * _lm_phi_face[j][qp];
324 : }
325 : }
326 839866 : }
327 :
328 : void
329 33439 : 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 100317 : for (const auto qp : make_range(qrule_face.n_points()))
339 : {
340 200634 : const auto scalar_value = dirichlet_value(
341 66878 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
342 66878 : _ti.determineState());
343 66878 : const auto qp_term = JxW_face[qp] * normals[qp] * scalar_value;
344 :
345 : // External boundary -> Dirichlet faces -> Vector equation RHS
346 434546 : for (const auto i : index_range(_qu_dof_indices))
347 367668 : vector_re(i) -= qp_term * _vector_phi_face[i][qp];
348 : }
349 33439 : }
350 :
351 : void
352 33439 : 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 100317 : for (const auto qp : make_range(qrule_face.n_points()))
364 : {
365 200634 : const auto scalar_value = dirichlet_value(
366 66878 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
367 66878 : _ti.determineState());
368 66878 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
369 66878 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
370 66878 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
371 66878 : const auto lm_qp_term = stab_qp_term * scalar_value;
372 :
373 224056 : for (const auto i : index_range(_u_dof_indices))
374 157178 : scalar_re(i) += (scalar_qp_term - vector_qp_term - lm_qp_term) * _scalar_phi_face[i][qp];
375 : }
376 33439 : }
377 :
378 : void
379 21966 : 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 65898 : for (const auto qp : make_range(qrule_face.n_points()))
386 : {
387 43932 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
388 43932 : const auto scalar_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
389 146272 : for (const auto i : index_range(_u_dof_indices))
390 : {
391 102340 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
392 777196 : for (const auto j : index_range(_qu_dof_indices))
393 674856 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
394 :
395 102340 : const auto scalar_qpi_term = scalar_qp_term * _scalar_phi_face[i][qp];
396 421904 : for (const auto j : index_range(_u_dof_indices))
397 319564 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
398 : }
399 : }
400 21966 : }
401 :
402 : void
403 33439 : 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 100317 : for (const auto qp : make_range(qrule.n_points()))
410 : {
411 66878 : const auto qp_term = JxW[qp] * sol[qp];
412 552370 : for (const auto i : index_range(phi))
413 485492 : re(i) -= phi[i][qp] * qp_term;
414 : }
415 33439 : }
416 :
417 : void
418 21966 : DiffusionLHDGAssemblyHelper::createIdentityJacobian(const MooseArray<Real> & JxW,
419 : const QBase & qrule,
420 : const MooseArray<std::vector<Real>> & phi,
421 : DenseMatrix<Number> & ke)
422 : {
423 65898 : for (const auto qp : make_range(qrule.n_points()))
424 363524 : for (const auto i : index_range(phi))
425 : {
426 319592 : const auto qpi_term = JxW[qp] * phi[i][qp];
427 2685144 : for (const auto j : index_range(phi))
428 2365552 : ke(i, j) -= phi[j][qp] * qpi_term;
429 : }
430 21966 : }
|
