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 43799 : DiffusionLHDGAssemblyHelper::validParams()
23 : {
24 43799 : auto params = emptyInputParameters();
25 43799 : params.addRequiredParam<NonlinearVariableName>(
26 : "gradient_variable", "The gradient of the diffusing specie concentration");
27 43799 : params.addRequiredParam<NonlinearVariableName>(
28 : "face_variable", "The concentration of the diffusing specie on faces");
29 43799 : params.addRequiredParam<MaterialPropertyName>("diffusivity", "The diffusivity");
30 131397 : params.addParam<Real>("tau",
31 87598 : 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 43799 : return params;
35 0 : }
36 :
37 504 : 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 504 : const THREAD_ID tid)
45 : : ADFunctorInterface(moose_obj),
46 504 : _u_var(sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("variable"))),
47 1008 : _grad_u_var(sys.getFieldVariable<RealVectorValue>(
48 504 : tid, moose_obj->getParam<NonlinearVariableName>("gradient_variable"))),
49 1008 : _u_face_var(sys.getFieldVariable<Real>(
50 504 : tid, moose_obj->getParam<NonlinearVariableName>("face_variable"))),
51 504 : _qu_dof_indices(_grad_u_var.dofIndices()),
52 504 : _u_dof_indices(_u_var.dofIndices()),
53 504 : _lm_u_dof_indices(_u_face_var.dofIndices()),
54 504 : _qu_sol(_grad_u_var.sln()),
55 504 : _u_sol(_u_var.sln()),
56 504 : _lm_u_sol(_u_face_var.sln()),
57 504 : _vector_phi(_grad_u_var.phi()),
58 504 : _scalar_phi(_u_var.phi()),
59 504 : _grad_scalar_phi(_u_var.gradPhi()),
60 504 : _div_vector_phi(_grad_u_var.divPhi()),
61 504 : _vector_phi_face(_grad_u_var.phiFace()),
62 504 : _scalar_phi_face(_u_var.phiFace()),
63 504 : _lm_phi_face(_u_face_var.phiFace()),
64 504 : _diff(mpi->getMaterialProperty<Real>("diffusivity")),
65 504 : _ti(*ti),
66 504 : _tau(moose_obj->getParam<Real>("tau")),
67 504 : _cached_elem(nullptr),
68 504 : _moose_obj(*moose_obj),
69 504 : _dhah_fe_problem(fe_problem),
70 1008 : _dhah_sys(sys)
71 : {
72 504 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<RealVectorValue> &>(_grad_u_var));
73 504 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_u_face_var));
74 504 : }
75 :
76 : void
77 504 : DiffusionLHDGAssemblyHelper::checkCoupling()
78 : {
79 : // This check must occur after FEProblemBase::init()
80 504 : if (_dhah_fe_problem.coupling() == Moose::COUPLING_FULL)
81 0 : return;
82 504 : else if (_dhah_fe_problem.coupling() == Moose::COUPLING_CUSTOM)
83 : {
84 504 : const auto * const cm = _dhah_fe_problem.couplingMatrix(_dhah_sys.number());
85 2016 : for (const auto i : make_range(cm->size()))
86 6048 : for (const auto j : make_range(cm->size()))
87 4536 : if ((*cm)(i, j) != true)
88 0 : goto error;
89 :
90 504 : 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 355368 : 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 1776840 : for (const auto qp : make_range(qrule.n_points()))
109 : {
110 1421472 : const auto vector_qp_term = JxW[qp] * vector_sol[qp];
111 1421472 : const auto scalar_qp_term = JxW[qp] * scalar_sol[qp];
112 9220512 : for (const auto i : index_range(vector_re))
113 : {
114 : // Vector equation dependence on vector dofs
115 7799040 : vector_re(i) += _vector_phi[i][qp] * vector_qp_term;
116 :
117 : // Vector equation dependence on scalar dofs
118 7799040 : vector_re(i) += _div_vector_phi[i][qp] * scalar_qp_term;
119 : }
120 : }
121 355368 : }
122 :
123 : void
124 235936 : DiffusionLHDGAssemblyHelper::vectorVolumeJacobian(const MooseArray<Real> & JxW,
125 : const QBase & qrule,
126 : DenseMatrix<Number> & vector_vector_jac,
127 : DenseMatrix<Number> & vector_scalar_jac)
128 : {
129 1179680 : for (const auto qp : make_range(qrule.n_points()))
130 6105568 : for (const auto i : make_range(vector_vector_jac.m()))
131 : {
132 : // Vector equation dependence on vector dofs
133 5161824 : const auto vector_qpi_term = JxW[qp] * _vector_phi[i][qp];
134 37076032 : for (const auto j : make_range(vector_vector_jac.n()))
135 31914208 : vector_vector_jac(i, j) += vector_qpi_term * _vector_phi[j][qp];
136 :
137 : // Vector equation dependence on scalar dofs
138 5161824 : const auto scalar_qpi_term = JxW[qp] * _div_vector_phi[i][qp];
139 19646016 : for (const auto j : make_range(vector_scalar_jac.n()))
140 14484192 : vector_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi[j][qp];
141 : }
142 235936 : }
143 :
144 : void
145 355368 : 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 1776840 : for (const auto qp : make_range(qrule.n_points()))
154 : {
155 1421472 : const auto vector_qp_term = JxW[qp] * _diff[qp] * vector_field[qp];
156 : // Evaluate source
157 : const auto f =
158 1421472 : source(Moose::ElemQpArg{current_elem, qp, &qrule, q_point[qp]}, _ti.determineState());
159 1421472 : const auto source_qp_term = JxW[qp] * f;
160 :
161 4740832 : for (const auto i : index_range(scalar_re))
162 : {
163 3319360 : scalar_re(i) += _grad_scalar_phi[i][qp] * vector_qp_term;
164 :
165 : // Scalar equation RHS
166 3319360 : scalar_re(i) -= _scalar_phi[i][qp] * source_qp_term;
167 : }
168 : }
169 355368 : }
170 :
171 : void
172 235936 : DiffusionLHDGAssemblyHelper::scalarVolumeJacobian(const MooseArray<Real> & JxW,
173 : const QBase & qrule,
174 : DenseMatrix<Number> & scalar_vector_jac)
175 : {
176 1179680 : for (const auto qp : make_range(qrule.n_points()))
177 : {
178 943744 : const auto qp_term = JxW[qp] * _diff[qp];
179 3135568 : for (const auto i : make_range(scalar_vector_jac.m()))
180 : {
181 2191824 : const auto qpi_term = qp_term * _grad_scalar_phi[i][qp];
182 : // Scalar equation dependence on vector dofs
183 16676016 : for (const auto j : make_range(scalar_vector_jac.n()))
184 14484192 : scalar_vector_jac(i, j) += qpi_term * _vector_phi[j][qp];
185 : }
186 : }
187 235936 : }
188 :
189 : void
190 1262104 : 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 3786312 : for (const auto qp : make_range(qrule_face.n_points()))
198 : {
199 2524208 : const auto qp_term = JxW_face[qp] * lm_sol[qp] * normals[qp];
200 16608368 : for (const auto i : index_range(vector_re))
201 14084160 : vector_re(i) -= _vector_phi_face[i][qp] * qp_term;
202 : }
203 1262104 : }
204 :
205 : void
206 838736 : 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 2516208 : for (const auto qp : make_range(qrule_face.n_points()))
212 : {
213 1677472 : const auto qp_term = JxW_face[qp] * normals[qp];
214 : // Vector equation dependence on LM dofs
215 11014528 : for (const auto i : make_range(vector_lm_jac.m()))
216 : {
217 9337056 : const auto qpi_term = qp_term * _vector_phi_face[i][qp];
218 79694496 : for (const auto j : make_range(vector_lm_jac.n()))
219 70357440 : vector_lm_jac(i, j) -= qpi_term * _lm_phi_face[j][qp];
220 : }
221 : }
222 838736 : }
223 :
224 : void
225 1262104 : 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 3786312 : for (const auto qp : make_range(qrule_face.n_points()))
234 : {
235 : // vector
236 2524208 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
237 : // stabilization term
238 2524208 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
239 : // scalar from stabilization term
240 2524208 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
241 : // lm from stabilization term
242 2524208 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
243 8493440 : for (const auto i : index_range(scalar_re))
244 5969232 : scalar_re(i) += _scalar_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
245 : }
246 1262104 : }
247 :
248 : void
249 838736 : 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 2516208 : for (const auto qp : make_range(qrule_face.n_points()))
257 : {
258 1677472 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
259 1677472 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
260 :
261 5627408 : for (const auto i : make_range(scalar_vector_jac.m()))
262 : {
263 3949936 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
264 30505552 : for (const auto j : make_range(scalar_vector_jac.n()))
265 26555616 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
266 :
267 3949936 : const auto scalar_qpi_term = stab_qp_term * _scalar_phi_face[i][qp];
268 16509152 : for (const auto j : make_range(scalar_scalar_jac.n()))
269 12559216 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
270 33622512 : for (const auto j : make_range(scalar_lm_jac.n()))
271 29672576 : scalar_lm_jac(i, j) -= scalar_qpi_term * _lm_phi_face[j][qp];
272 : }
273 : }
274 838736 : }
275 :
276 : void
277 1262104 : 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 3786312 : for (const auto qp : make_range(qrule_face.n_points()))
286 : {
287 : // vector
288 2524208 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
289 : // stabilization term
290 2524208 : const auto stab_qp_term = JxW_face[qp] * _tau * (normals[qp] * normals[qp]);
291 : // scalar from stabilization term
292 2524208 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
293 : // lm from stabilization term
294 2524208 : const auto lm_qp_term = stab_qp_term * lm_sol[qp];
295 21253456 : for (const auto i : index_range(lm_re))
296 18729248 : lm_re(i) += _lm_phi_face[i][qp] * (scalar_qp_term - vector_qp_term - lm_qp_term);
297 : }
298 1262104 : }
299 :
300 : void
301 838736 : 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 2516208 : for (const auto qp : make_range(qrule_face.n_points()))
309 : {
310 1677472 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
311 1677472 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
312 :
313 14131488 : for (const auto i : make_range(lm_vec_jac.m()))
314 : {
315 12454016 : const auto vector_qpi_term = vector_qp_term * _lm_phi_face[i][qp];
316 82811456 : for (const auto j : make_range(lm_vec_jac.n()))
317 70357440 : lm_vec_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
318 :
319 12454016 : const auto lm_qpi_term = stab_qp_term * _lm_phi_face[i][qp];
320 42126592 : for (const auto j : make_range(lm_scalar_jac.n()))
321 29672576 : lm_scalar_jac(i, j) += lm_qpi_term * _scalar_phi_face[j][qp];
322 106291584 : for (const auto j : make_range(lm_lm_jac.n()))
323 93837568 : lm_lm_jac(i, j) -= lm_qpi_term * _lm_phi_face[j][qp];
324 : }
325 : }
326 838736 : }
327 :
328 : void
329 33264 : 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 99792 : for (const auto qp : make_range(qrule_face.n_points()))
339 : {
340 199584 : const auto scalar_value = dirichlet_value(
341 66528 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
342 66528 : _ti.determineState());
343 66528 : const auto qp_term = JxW_face[qp] * normals[qp] * scalar_value;
344 :
345 : // External boundary -> Dirichlet faces -> Vector equation RHS
346 432096 : for (const auto i : index_range(_qu_dof_indices))
347 365568 : vector_re(i) -= qp_term * _vector_phi_face[i][qp];
348 : }
349 33264 : }
350 :
351 : void
352 33264 : 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 99792 : for (const auto qp : make_range(qrule_face.n_points()))
364 : {
365 199584 : const auto scalar_value = dirichlet_value(
366 66528 : Moose::ElemSideQpArg{current_elem, current_side, qp, &qrule_face, q_point_face[qp]},
367 66528 : _ti.determineState());
368 66528 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * (vector_sol[qp] * normals[qp]);
369 66528 : const auto stab_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
370 66528 : const auto scalar_qp_term = stab_qp_term * scalar_sol[qp];
371 66528 : const auto lm_qp_term = stab_qp_term * scalar_value;
372 :
373 222656 : for (const auto i : index_range(_u_dof_indices))
374 156128 : scalar_re(i) += (scalar_qp_term - vector_qp_term - lm_qp_term) * _scalar_phi_face[i][qp];
375 : }
376 33264 : }
377 :
378 : void
379 21896 : 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 65688 : for (const auto qp : make_range(qrule_face.n_points()))
386 : {
387 43792 : const auto vector_qp_term = JxW_face[qp] * _diff[qp] * normals[qp];
388 43792 : const auto scalar_qp_term = JxW_face[qp] * _tau * normals[qp] * normals[qp];
389 145712 : for (const auto i : index_range(_u_dof_indices))
390 : {
391 101920 : const auto vector_qpi_term = vector_qp_term * _scalar_phi_face[i][qp];
392 774256 : for (const auto j : index_range(_qu_dof_indices))
393 672336 : scalar_vector_jac(i, j) -= vector_qpi_term * _vector_phi_face[j][qp];
394 :
395 101920 : const auto scalar_qpi_term = scalar_qp_term * _scalar_phi_face[i][qp];
396 420224 : for (const auto j : index_range(_u_dof_indices))
397 318304 : scalar_scalar_jac(i, j) += scalar_qpi_term * _scalar_phi_face[j][qp];
398 : }
399 : }
400 21896 : }
401 :
402 : void
403 33264 : 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 99792 : for (const auto qp : make_range(qrule.n_points()))
410 : {
411 66528 : const auto qp_term = JxW[qp] * sol[qp];
412 549920 : for (const auto i : index_range(phi))
413 483392 : re(i) -= phi[i][qp] * qp_term;
414 : }
415 33264 : }
416 :
417 : void
418 21896 : DiffusionLHDGAssemblyHelper::createIdentityJacobian(const MooseArray<Real> & JxW,
419 : const QBase & qrule,
420 : const MooseArray<std::vector<Real>> & phi,
421 : DenseMatrix<Number> & ke)
422 : {
423 65688 : for (const auto qp : make_range(qrule.n_points()))
424 362544 : for (const auto i : index_range(phi))
425 : {
426 318752 : const auto qpi_term = JxW[qp] * phi[i][qp];
427 2679264 : for (const auto j : index_range(phi))
428 2360512 : ke(i, j) -= phi[j][qp] * qpi_term;
429 : }
430 21896 : }
|
