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 "NavierStokesLHDGAssemblyHelper.h"
11 : #include "MooseVariableFE.h"
12 : #include "MooseVariableScalar.h"
13 : #include "SystemBase.h"
14 : #include "MooseMesh.h"
15 : #include "MooseObject.h"
16 : #include "NS.h"
17 : #include "TransientInterface.h"
18 :
19 : using namespace libMesh;
20 :
21 : InputParameters
22 1664 : NavierStokesLHDGAssemblyHelper::validParams()
23 : {
24 1664 : auto params = DiffusionLHDGAssemblyHelper::validParams();
25 1664 : params.addRequiredParam<NonlinearVariableName>(NS::pressure, "The pressure variable.");
26 3328 : params.addRequiredParam<NonlinearVariableName>("v", "The y-component of velocity");
27 3328 : params.addParam<NonlinearVariableName>("w", "The z-component of velocity");
28 3328 : params.renameParam("gradient_variable", "grad_u", "The gradient of the x-component of velocity");
29 3328 : params.addRequiredParam<NonlinearVariableName>("grad_v",
30 : "The gradient of the y-component of velocity");
31 3328 : params.addParam<NonlinearVariableName>("grad_w", "The gradient of the z-component of velocity");
32 3328 : params.renameParam("face_variable", "face_u", "The x-component of the face velocity");
33 3328 : params.addRequiredParam<NonlinearVariableName>("face_v", "The y-component of the face velocity");
34 3328 : params.addParam<NonlinearVariableName>("face_w", "The z-component of the face velocity");
35 3328 : params.addParam<NonlinearVariableName>(
36 : "enclosure_lm",
37 : "For enclosed problems like the lid driven cavity this variable can be provided to remove "
38 : "the pressure nullspace");
39 3328 : params.renameParam("diffusivity", NS::mu, "The dynamic viscosity");
40 1664 : params.addRequiredParam<Real>(NS::density, "The density");
41 1664 : return params;
42 0 : }
43 :
44 832 : NavierStokesLHDGAssemblyHelper::NavierStokesLHDGAssemblyHelper(
45 : const MooseObject * const moose_obj,
46 : MaterialPropertyInterface * const mpi,
47 : MooseVariableDependencyInterface * const mvdi,
48 : const TransientInterface * const ti,
49 : const FEProblemBase & fe_problem,
50 : SystemBase & sys,
51 : const MooseMesh & mesh,
52 832 : const THREAD_ID tid)
53 : : DiffusionLHDGAssemblyHelper(moose_obj, mpi, mvdi, ti, fe_problem, sys, tid),
54 : // vars
55 832 : _v_var(sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("v"))),
56 832 : _w_var(mesh.dimension() > 2
57 832 : ? &sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("w"))
58 : : nullptr),
59 832 : _grad_v_var(sys.getFieldVariable<RealVectorValue>(
60 832 : tid, moose_obj->getParam<NonlinearVariableName>("grad_v"))),
61 832 : _grad_w_var(mesh.dimension() > 2
62 832 : ? &sys.getFieldVariable<RealVectorValue>(
63 832 : tid, moose_obj->getParam<NonlinearVariableName>("grad_w"))
64 : : nullptr),
65 832 : _v_face_var(
66 1664 : sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("face_v"))),
67 832 : _w_face_var(
68 832 : mesh.dimension() > 2
69 832 : ? &sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("face_w"))
70 : : nullptr),
71 832 : _pressure_var(
72 832 : sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>(NS::pressure))),
73 832 : _enclosure_lm_var(moose_obj->isParamValid("enclosure_lm")
74 1168 : ? &sys.getScalarVariable(
75 1504 : tid, moose_obj->getParam<NonlinearVariableName>("enclosure_lm"))
76 : : nullptr),
77 : // dof indices
78 832 : _qv_dof_indices(_grad_v_var.dofIndices()),
79 832 : _v_dof_indices(_v_var.dofIndices()),
80 832 : _lm_v_dof_indices(_v_face_var.dofIndices()),
81 832 : _qw_dof_indices(_grad_w_var ? &_grad_w_var->dofIndices() : nullptr),
82 832 : _w_dof_indices(_w_var ? &_w_var->dofIndices() : nullptr),
83 832 : _lm_w_dof_indices(_w_face_var ? &_w_face_var->dofIndices() : nullptr),
84 832 : _p_dof_indices(_pressure_var.dofIndices()),
85 832 : _global_lm_dof_indices(_enclosure_lm_var ? &_enclosure_lm_var->dofIndices() : nullptr),
86 : // solutions
87 832 : _qv_sol(_grad_v_var.sln()),
88 832 : _v_sol(_v_var.sln()),
89 832 : _lm_v_sol(_v_face_var.sln()),
90 832 : _qw_sol(_grad_w_var ? &_grad_w_var->sln() : nullptr),
91 832 : _w_sol(_w_var ? &_w_var->sln() : nullptr),
92 832 : _lm_w_sol(_w_face_var ? &_w_face_var->sln() : nullptr),
93 832 : _p_sol(_pressure_var.sln()),
94 1168 : _global_lm_dof_value(_enclosure_lm_var ? &_enclosure_lm_var->sln() : nullptr),
95 1664 : _rho(moose_obj->getParam<Real>(NS::density))
96 : {
97 832 : if (mesh.dimension() > 2)
98 0 : mooseError("3D not yet implemented");
99 :
100 832 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_v_var));
101 832 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<RealVectorValue> &>(_grad_v_var));
102 832 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_v_face_var));
103 832 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_pressure_var));
104 :
105 832 : const auto vel_type = _u_var.feType();
106 1664 : auto check_type = [&vel_type](const auto & var)
107 : {
108 1664 : if (vel_type != var.feType())
109 0 : mooseError(
110 : var.name(),
111 : "does not have the same finite element type as the x-component velocity. All scalar "
112 : "field finite element types must be the same for the Navier-Stokes L-HDG implementation");
113 1664 : };
114 832 : check_type(_v_var);
115 832 : if (_w_var)
116 0 : check_type(*_w_var);
117 832 : check_type(_pressure_var);
118 :
119 832 : const auto grad_type = _grad_u_var.feType();
120 832 : auto check_grad_type = [&grad_type](const auto & var)
121 : {
122 832 : if (grad_type != var.feType())
123 0 : mooseError(
124 : var.name(),
125 : "does not have the same finite element type as the x-component velocity gradient. All "
126 : "vector field finite element types must be the same for the Navier-Stokes L-HDG "
127 : "implementation");
128 832 : };
129 832 : check_grad_type(_grad_v_var);
130 832 : if (_grad_w_var)
131 0 : check_grad_type(*_grad_w_var);
132 :
133 832 : const auto lm_type = _u_face_var.feType();
134 832 : auto check_lm_type = [&lm_type](const auto & var)
135 : {
136 832 : if (lm_type != var.feType())
137 0 : mooseError(var.name(),
138 : "does not have the same finite element type as the x-component face variable. All "
139 : "face variable finite element types must be the same for the Navier-Stokes L-HDG "
140 : "implementation");
141 832 : };
142 832 : check_lm_type(_v_face_var);
143 832 : if (_w_face_var)
144 0 : check_lm_type(*_w_face_var);
145 832 : }
146 :
147 : void
148 3460992 : NavierStokesLHDGAssemblyHelper::scalarVolumeResidual(const MooseArray<Gradient> & vel_gradient,
149 : const unsigned int vel_component,
150 : const Moose::Functor<Real> & body_force,
151 : const MooseArray<Real> & JxW,
152 : const QBase & qrule,
153 : const Elem * const current_elem,
154 : const MooseArray<Point> & q_point,
155 : DenseVector<Number> & scalar_re)
156 : {
157 3460992 : DiffusionLHDGAssemblyHelper::scalarVolumeResidual(
158 : vel_gradient, body_force, JxW, qrule, current_elem, q_point, scalar_re);
159 :
160 17304960 : for (const auto qp : make_range(qrule.n_points()))
161 : {
162 13843968 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_u_sol, _v_sol, qp, vel_component);
163 : Gradient qp_p;
164 13843968 : qp_p(vel_component) = _p_sol[qp];
165 :
166 55375872 : for (const auto i : index_range(scalar_re))
167 : {
168 : // Scalar equation dependence on pressure dofs
169 41531904 : scalar_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * qp_p);
170 :
171 : // Scalar equation dependence on scalar dofs
172 41531904 : scalar_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * rho_vel_cross_vel);
173 : }
174 : }
175 3460992 : }
176 :
177 : void
178 985712 : NavierStokesLHDGAssemblyHelper::scalarVolumeJacobian(const unsigned int vel_component,
179 : const MooseArray<Real> & JxW,
180 : const QBase & qrule,
181 : DenseMatrix<Number> & scalar_vector_jac,
182 : DenseMatrix<Number> & scalar_u_vel_jac,
183 : DenseMatrix<Number> & scalar_v_vel_jac,
184 : DenseMatrix<Number> & scalar_p_jac)
185 : {
186 985712 : DiffusionLHDGAssemblyHelper::scalarVolumeJacobian(JxW, qrule, scalar_vector_jac);
187 :
188 3942848 : for (const auto i : make_range(scalar_vector_jac.m()))
189 14785680 : for (const auto qp : make_range(qrule.n_points()))
190 : {
191 : // Scalar equation dependence on pressure dofs
192 47314176 : for (const auto j : make_range(scalar_p_jac.n()))
193 : {
194 : Gradient p_phi;
195 35485632 : p_phi(vel_component) = _scalar_phi[j][qp];
196 35485632 : scalar_p_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * p_phi);
197 : }
198 :
199 : // Scalar equation dependence on scalar dofs
200 : mooseAssert(scalar_u_vel_jac.n() == scalar_v_vel_jac.n(), "These must be the same size");
201 47314176 : for (const auto j : make_range(scalar_u_vel_jac.n()))
202 : {
203 : // derivatives wrt 0th component of velocity
204 : {
205 : const auto rho_vel_cross_vel =
206 35485632 : rhoVelCrossVelJacobian(_u_sol, _v_sol, qp, vel_component, 0, _scalar_phi, j);
207 35485632 : scalar_u_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * rho_vel_cross_vel);
208 : }
209 : // derivatives wrt 1th component of velocity
210 : {
211 : const auto rho_vel_cross_vel =
212 35485632 : rhoVelCrossVelJacobian(_u_sol, _v_sol, qp, vel_component, 1, _scalar_phi, j);
213 35485632 : scalar_v_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * rho_vel_cross_vel);
214 : }
215 : }
216 : }
217 985712 : }
218 :
219 : void
220 1730496 : NavierStokesLHDGAssemblyHelper::pressureVolumeResidual(
221 : const Moose::Functor<Real> & pressure_mms_forcing_function,
222 : const MooseArray<Real> & JxW,
223 : const QBase & qrule,
224 : const Elem * const current_elem,
225 : const MooseArray<Point> & q_point,
226 : DenseVector<Number> & pressure_re,
227 : DenseVector<Number> & global_lm_re)
228 : {
229 8652480 : for (const auto qp : make_range(qrule.n_points()))
230 : {
231 : // Prepare forcing function
232 6921984 : const auto f = pressure_mms_forcing_function(
233 6921984 : Moose::ElemQpArg{current_elem, qp, &qrule, q_point[qp]}, _ti.determineState());
234 :
235 6921984 : const Gradient vel(_u_sol[qp], _v_sol[qp]);
236 27687936 : for (const auto i : make_range(pressure_re.size()))
237 : {
238 20765952 : pressure_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * vel);
239 :
240 : // Pressure equation forcing function RHS
241 20765952 : pressure_re(i) -= JxW[qp] * _scalar_phi[i][qp] * f;
242 :
243 20765952 : if (_enclosure_lm_var)
244 : {
245 : mooseAssert(
246 : _global_lm_dof_value->size() == 1,
247 : "There should only be one degree of freedom for removing the pressure nullspace");
248 9973824 : pressure_re(i) -= JxW[qp] * _scalar_phi[i][qp] * (*_global_lm_dof_value)[0];
249 : }
250 : }
251 :
252 6921984 : if (_enclosure_lm_var)
253 3324608 : global_lm_re(0) -= JxW[qp] * _p_sol[qp];
254 : }
255 1730496 : }
256 :
257 : void
258 492856 : NavierStokesLHDGAssemblyHelper::pressureVolumeJacobian(const MooseArray<Real> & JxW,
259 : const QBase & qrule,
260 : DenseMatrix<Number> & p_u_vel_jac,
261 : DenseMatrix<Number> & p_v_vel_jac,
262 : DenseMatrix<Number> & p_global_lm_jac,
263 : DenseMatrix<Number> & global_lm_p_jac)
264 : {
265 2464280 : for (const auto qp : make_range(qrule.n_points()))
266 : {
267 7885696 : for (const auto i : make_range(p_u_vel_jac.m()))
268 : {
269 23657088 : for (const auto j : make_range(p_u_vel_jac.n()))
270 : {
271 : {
272 17742816 : const Gradient phi(_scalar_phi[j][qp], 0);
273 17742816 : p_u_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * phi);
274 : }
275 : {
276 : const Gradient phi(0, _scalar_phi[j][qp]);
277 17742816 : p_v_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * phi);
278 : }
279 : }
280 5914272 : if (_enclosure_lm_var)
281 : p_global_lm_jac(i, 0) -= JxW[qp] * _scalar_phi[i][qp];
282 : }
283 :
284 1971424 : if (_enclosure_lm_var)
285 : {
286 3478528 : for (const auto j : make_range(global_lm_p_jac.n()))
287 2608896 : global_lm_p_jac(0, j) -= JxW[qp] * _scalar_phi[j][qp];
288 : }
289 : }
290 492856 : }
291 :
292 : RealVectorValue
293 54731552 : NavierStokesLHDGAssemblyHelper::rhoVelCrossVelResidual(const MooseArray<Number> & u_sol,
294 : const MooseArray<Number> & v_sol,
295 : const unsigned int qp,
296 : const unsigned int vel_component)
297 : {
298 : const RealVectorValue U(u_sol[qp], v_sol[qp]);
299 54731552 : return _rho * U * U(vel_component);
300 : }
301 :
302 : RealVectorValue
303 698851296 : NavierStokesLHDGAssemblyHelper::rhoVelCrossVelJacobian(const MooseArray<Number> & u_sol,
304 : const MooseArray<Number> & v_sol,
305 : const unsigned int qp,
306 : const unsigned int vel_component,
307 : const unsigned int vel_j_component,
308 : const MooseArray<std::vector<Real>> & phi,
309 : const unsigned int j)
310 : {
311 : const RealVectorValue U(u_sol[qp], v_sol[qp]);
312 : RealVectorValue vector_phi;
313 698851296 : vector_phi(vel_j_component) = phi[j][qp];
314 : auto ret = vector_phi * U(vel_component);
315 698851296 : if (vel_component == vel_j_component)
316 : ret += U * phi[j][qp];
317 : ret *= _rho;
318 698851296 : return ret;
319 : }
320 :
321 : void
322 5110948 : NavierStokesLHDGAssemblyHelper::pressureFaceResidual(const MooseArray<Real> & JxW_face,
323 : const QBase & qrule_face,
324 : const MooseArray<Point> & normals,
325 : DenseVector<Number> & pressure_re)
326 : {
327 15332844 : for (const auto qp : make_range(qrule_face.n_points()))
328 : {
329 10221896 : const Gradient vel(_lm_u_sol[qp], _lm_v_sol[qp]);
330 : const auto vdotn = vel * normals[qp];
331 40887584 : for (const auto i : make_range(pressure_re.size()))
332 30665688 : pressure_re(i) += JxW_face[qp] * vdotn * _scalar_phi_face[i][qp];
333 : }
334 5110948 : }
335 :
336 : void
337 1455574 : NavierStokesLHDGAssemblyHelper::pressureFaceJacobian(const MooseArray<Real> & JxW_face,
338 : const QBase & qrule_face,
339 : const MooseArray<Point> & normals,
340 : DenseMatrix<Number> & p_lm_u_vel_jac,
341 : DenseMatrix<Number> & p_lm_v_vel_jac)
342 : {
343 : mooseAssert((p_lm_u_vel_jac.m() == p_lm_v_vel_jac.m()) &&
344 : (p_lm_u_vel_jac.n() == p_lm_v_vel_jac.n()),
345 : "We already checked that lm finite element types are the same, so these matrices "
346 : "should be the same size");
347 5822296 : for (const auto i : make_range(p_lm_u_vel_jac.m()))
348 30567054 : for (const auto j : make_range(p_lm_u_vel_jac.n()))
349 78600996 : for (const auto qp : make_range(qrule_face.n_points()))
350 : {
351 : {
352 52400664 : const Gradient phi(_lm_phi_face[j][qp], 0);
353 52400664 : p_lm_u_vel_jac(i, j) += JxW_face[qp] * phi * normals[qp] * _scalar_phi_face[i][qp];
354 : }
355 : {
356 : const Gradient phi(0, _lm_phi_face[j][qp]);
357 52400664 : p_lm_v_vel_jac(i, j) += JxW_face[qp] * phi * normals[qp] * _scalar_phi_face[i][qp];
358 : }
359 : }
360 1455574 : }
361 :
362 : void
363 10221896 : NavierStokesLHDGAssemblyHelper::scalarFaceResidual(const MooseArray<Gradient> & vector_sol,
364 : const MooseArray<Number> & scalar_sol,
365 : const MooseArray<Number> & lm_sol,
366 : const unsigned int vel_component,
367 : const MooseArray<Real> & JxW_face,
368 : const QBase & qrule_face,
369 : const MooseArray<Point> & normals,
370 : DenseVector<Number> & scalar_re)
371 : {
372 10221896 : DiffusionLHDGAssemblyHelper::scalarFaceResidual(
373 : vector_sol, scalar_sol, lm_sol, JxW_face, qrule_face, normals, scalar_re);
374 :
375 30665688 : for (const auto qp : make_range(qrule_face.n_points()))
376 : {
377 : Gradient qp_p;
378 20443792 : qp_p(vel_component) = _p_sol[qp];
379 20443792 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_lm_u_sol, _lm_v_sol, qp, vel_component);
380 :
381 81775168 : for (const auto i : make_range(scalar_re.size()))
382 : {
383 : // pressure
384 61331376 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * (qp_p * normals[qp]);
385 :
386 : // lm from convection term
387 61331376 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
388 : }
389 : }
390 10221896 : }
391 :
392 : void
393 2911148 : NavierStokesLHDGAssemblyHelper::scalarFaceJacobian(const unsigned int vel_component,
394 : const MooseArray<Real> & JxW_face,
395 : const QBase & qrule_face,
396 : const MooseArray<Point> & normals,
397 : DenseMatrix<Number> & scalar_vector_jac,
398 : DenseMatrix<Number> & scalar_scalar_jac,
399 : DenseMatrix<Number> & scalar_lm_jac,
400 : DenseMatrix<Number> & scalar_p_jac,
401 : DenseMatrix<Number> & scalar_lm_u_vel_jac,
402 : DenseMatrix<Number> & scalar_lm_v_vel_jac)
403 :
404 : {
405 2911148 : DiffusionLHDGAssemblyHelper::scalarFaceJacobian(
406 : JxW_face, qrule_face, normals, scalar_vector_jac, scalar_scalar_jac, scalar_lm_jac);
407 :
408 11644592 : for (const auto i : make_range(scalar_lm_u_vel_jac.m()))
409 26200332 : for (const auto qp : make_range(qrule_face.n_points()))
410 : {
411 69867552 : for (const auto j : make_range(scalar_p_jac.n()))
412 : {
413 : Gradient p_phi;
414 52400664 : p_phi(vel_component) = _scalar_phi_face[j][qp];
415 : // pressure
416 52400664 : scalar_p_jac(i, j) += JxW_face[qp] * _scalar_phi_face[i][qp] * (p_phi * normals[qp]);
417 : }
418 :
419 122268216 : for (const auto j : make_range(scalar_lm_u_vel_jac.n()))
420 : {
421 : //
422 : // from convection term
423 : //
424 :
425 : // derivatives wrt 0th component of velocity
426 : {
427 : const auto rho_vel_cross_vel =
428 104801328 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 0, _lm_phi_face, j);
429 104801328 : scalar_lm_u_vel_jac(i, j) +=
430 104801328 : JxW_face[qp] * _scalar_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
431 : }
432 : // derivatives wrt 1th component of velocity
433 : {
434 : const auto rho_vel_cross_vel =
435 104801328 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 1, _lm_phi_face, j);
436 104801328 : scalar_lm_v_vel_jac(i, j) +=
437 104801328 : JxW_face[qp] * _scalar_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
438 : }
439 : }
440 : }
441 2911148 : }
442 :
443 : void
444 10221896 : NavierStokesLHDGAssemblyHelper::lmFaceResidual(const MooseArray<Gradient> & vector_sol,
445 : const MooseArray<Number> & scalar_sol,
446 : const MooseArray<Number> & lm_sol,
447 : const unsigned int vel_component,
448 : const MooseArray<Real> & JxW_face,
449 : const QBase & qrule_face,
450 : const MooseArray<Point> & normals,
451 : const Elem * const neigh,
452 : DenseVector<Number> & lm_re)
453 : {
454 10221896 : DiffusionLHDGAssemblyHelper::lmFaceResidual(
455 : vector_sol, scalar_sol, lm_sol, JxW_face, qrule_face, normals, lm_re);
456 :
457 30665688 : for (const auto qp : make_range(qrule_face.n_points()))
458 : {
459 : Gradient qp_p;
460 20443792 : qp_p(vel_component) = _p_sol[qp];
461 20443792 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_lm_u_sol, _lm_v_sol, qp, vel_component);
462 :
463 143106544 : for (const auto i : make_range(lm_re.size()))
464 : {
465 : // pressure
466 122662752 : lm_re(i) += JxW_face[qp] * _lm_phi_face[i][qp] * (qp_p * normals[qp]);
467 :
468 : // If we are an internal face we add the convective term. On the outflow boundary we do not
469 : // zero out the convection term, e.g. we are going to set q + p + tau * (u - u_hat) to zero
470 122662752 : if (neigh)
471 : // lm from convection term
472 122409984 : lm_re(i) += JxW_face[qp] * _lm_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
473 : }
474 : }
475 10221896 : }
476 :
477 : void
478 2911148 : NavierStokesLHDGAssemblyHelper::lmFaceJacobian(const unsigned int vel_component,
479 : const MooseArray<Real> & JxW_face,
480 : const QBase & qrule_face,
481 : const MooseArray<Point> & normals,
482 : const Elem * const neigh,
483 : DenseMatrix<Number> & lm_vec_jac,
484 : DenseMatrix<Number> & lm_scalar_jac,
485 : DenseMatrix<Number> & lm_lm_jac,
486 : DenseMatrix<Number> & lm_p_jac,
487 : DenseMatrix<Number> & lm_lm_u_vel_jac,
488 : DenseMatrix<Number> & lm_lm_v_vel_jac)
489 : {
490 2911148 : DiffusionLHDGAssemblyHelper::lmFaceJacobian(
491 : JxW_face, qrule_face, normals, lm_vec_jac, lm_scalar_jac, lm_lm_jac);
492 :
493 20378036 : for (const auto i : make_range(lm_p_jac.m()))
494 52400664 : for (const auto qp : make_range(qrule_face.n_points()))
495 : {
496 139735104 : for (const auto j : make_range(lm_p_jac.n()))
497 : {
498 : Gradient p_phi;
499 104801328 : p_phi(vel_component) = _scalar_phi_face[j][qp];
500 104801328 : lm_p_jac(i, j) += JxW_face[qp] * _lm_phi_face[i][qp] * (p_phi * normals[qp]);
501 : }
502 :
503 244536432 : for (const auto j : make_range(lm_lm_u_vel_jac.n()))
504 209602656 : if (neigh)
505 : {
506 : // derivatives wrt 0th component of velocity
507 : {
508 : const auto rho_vel_cross_vel =
509 209138688 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 0, _lm_phi_face, j);
510 209138688 : lm_lm_u_vel_jac(i, j) +=
511 209138688 : JxW_face[qp] * _lm_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
512 : }
513 : // derivatives wrt 1th component of velocity
514 : {
515 : const auto rho_vel_cross_vel =
516 209138688 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 1, _lm_phi_face, j);
517 209138688 : lm_lm_v_vel_jac(i, j) +=
518 209138688 : JxW_face[qp] * _lm_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
519 : }
520 : }
521 : }
522 2911148 : }
523 :
524 : void
525 80540 : NavierStokesLHDGAssemblyHelper::pressureDirichletResidual(
526 : const std::array<const Moose::Functor<Real> *, 3> & dirichlet_vel,
527 : const MooseArray<Real> & JxW_face,
528 : const QBase & qrule_face,
529 : const MooseArray<Point> & normals,
530 : const Elem * const current_elem,
531 : const unsigned int current_side,
532 : const MooseArray<Point> & q_point_face,
533 : DenseVector<Number> & pressure_re)
534 : {
535 241620 : for (const auto qp : make_range(qrule_face.n_points()))
536 : {
537 : const Moose::ElemSideQpArg elem_side_qp_arg{
538 161080 : current_elem, current_side, qp, &qrule_face, q_point_face[qp]};
539 161080 : const auto time_arg = _ti.determineState();
540 161080 : const RealVectorValue dirichlet_velocity((*dirichlet_vel[0])(elem_side_qp_arg, time_arg),
541 161080 : (*dirichlet_vel[1])(elem_side_qp_arg, time_arg),
542 161080 : (*dirichlet_vel[2])(elem_side_qp_arg, time_arg));
543 : const auto vdotn = dirichlet_velocity * normals[qp];
544 644320 : for (const auto i : make_range(pressure_re.size()))
545 483240 : pressure_re(i) += JxW_face[qp] * vdotn * _scalar_phi_face[i][qp];
546 : }
547 80540 : }
548 :
549 : void
550 161080 : NavierStokesLHDGAssemblyHelper::scalarDirichletResidual(
551 : const MooseArray<Gradient> & vector_sol,
552 : const MooseArray<Number> & scalar_sol,
553 : const unsigned int vel_component,
554 : const std::array<const Moose::Functor<Real> *, 3> & dirichlet_vel,
555 : const MooseArray<Real> & JxW_face,
556 : const QBase & qrule_face,
557 : const MooseArray<Point> & normals,
558 : const Elem * const current_elem,
559 : const unsigned int current_side,
560 : const MooseArray<Point> & q_point_face,
561 : DenseVector<Number> & scalar_re)
562 : {
563 161080 : DiffusionLHDGAssemblyHelper::scalarDirichletResidual(vector_sol,
564 : scalar_sol,
565 161080 : *dirichlet_vel[vel_component],
566 : JxW_face,
567 : qrule_face,
568 : normals,
569 : current_elem,
570 : current_side,
571 : q_point_face,
572 : scalar_re);
573 :
574 483240 : for (const auto qp : make_range(qrule_face.n_points()))
575 : {
576 : Gradient qp_p;
577 322160 : qp_p(vel_component) = _p_sol[qp];
578 :
579 : const Moose::ElemSideQpArg elem_side_qp_arg{
580 322160 : current_elem, current_side, qp, &qrule_face, q_point_face[qp]};
581 322160 : const auto time_arg = _ti.determineState();
582 322160 : const RealVectorValue dirichlet_velocity((*dirichlet_vel[0])(elem_side_qp_arg, time_arg),
583 322160 : (*dirichlet_vel[1])(elem_side_qp_arg, time_arg),
584 322160 : (*dirichlet_vel[2])(elem_side_qp_arg, time_arg));
585 322160 : const auto scalar_value = dirichlet_velocity(vel_component);
586 :
587 1288640 : for (const auto i : make_range(scalar_re.size()))
588 : {
589 : // pressure
590 966480 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * (qp_p * normals[qp]);
591 :
592 : // dirichlet lm from advection term
593 966480 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] *
594 966480 : (_rho * dirichlet_velocity * normals[qp]) * scalar_value;
595 : }
596 : }
597 161080 : }
598 :
599 : void
600 45988 : NavierStokesLHDGAssemblyHelper::scalarDirichletJacobian(const unsigned int vel_component,
601 : const MooseArray<Real> & JxW_face,
602 : const QBase & qrule_face,
603 : const MooseArray<Point> & normals,
604 : DenseMatrix<Number> & scalar_vector_jac,
605 : DenseMatrix<Number> & scalar_scalar_jac,
606 : DenseMatrix<Number> & scalar_pressure_jac)
607 : {
608 45988 : DiffusionLHDGAssemblyHelper::scalarDirichletJacobian(
609 : JxW_face, qrule_face, normals, scalar_vector_jac, scalar_scalar_jac);
610 :
611 183952 : for (const auto i : make_range(scalar_pressure_jac.m()))
612 551856 : for (const auto j : make_range(scalar_pressure_jac.n()))
613 1241676 : for (const auto qp : make_range(qrule_face.n_points()))
614 : {
615 : Gradient p_phi;
616 827784 : p_phi(vel_component) = _scalar_phi_face[j][qp];
617 : // pressure
618 827784 : scalar_pressure_jac(i, j) += JxW_face[qp] * _scalar_phi_face[i][qp] * (p_phi * normals[qp]);
619 : }
620 45988 : }
|
