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 2932 : NavierStokesLHDGAssemblyHelper::validParams()
23 : {
24 2932 : auto params = DiffusionLHDGAssemblyHelper::validParams();
25 2932 : params.addRequiredParam<NonlinearVariableName>(NS::pressure, "The pressure variable.");
26 5864 : params.addRequiredParam<NonlinearVariableName>("v", "The y-component of velocity");
27 5864 : params.addParam<NonlinearVariableName>("w", "The z-component of velocity");
28 5864 : params.renameParam("gradient_variable", "grad_u", "The gradient of the x-component of velocity");
29 5864 : params.addRequiredParam<NonlinearVariableName>("grad_v",
30 : "The gradient of the y-component of velocity");
31 5864 : params.addParam<NonlinearVariableName>("grad_w", "The gradient of the z-component of velocity");
32 5864 : params.renameParam("face_variable", "face_u", "The x-component of the face velocity");
33 5864 : params.addRequiredParam<NonlinearVariableName>("face_v", "The y-component of the face velocity");
34 5864 : params.addParam<NonlinearVariableName>("face_w", "The z-component of the face velocity");
35 5864 : 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 5864 : params.renameParam("diffusivity", NS::mu, "The dynamic viscosity");
40 2932 : params.addRequiredParam<Real>(NS::density, "The density");
41 2932 : return params;
42 0 : }
43 :
44 1466 : 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 1466 : const THREAD_ID tid)
53 : : DiffusionLHDGAssemblyHelper(moose_obj, mpi, mvdi, ti, fe_problem, sys, tid),
54 : // vars
55 1466 : _v_var(sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("v"))),
56 1466 : _w_var(mesh.dimension() > 2
57 1466 : ? &sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("w"))
58 : : nullptr),
59 1466 : _grad_v_var(sys.getFieldVariable<RealVectorValue>(
60 1466 : tid, moose_obj->getParam<NonlinearVariableName>("grad_v"))),
61 1466 : _grad_w_var(mesh.dimension() > 2
62 1466 : ? &sys.getFieldVariable<RealVectorValue>(
63 1466 : tid, moose_obj->getParam<NonlinearVariableName>("grad_w"))
64 : : nullptr),
65 1466 : _v_face_var(
66 2932 : sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("face_v"))),
67 1466 : _w_face_var(
68 1466 : mesh.dimension() > 2
69 1466 : ? &sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>("face_w"))
70 : : nullptr),
71 1466 : _pressure_var(
72 1466 : sys.getFieldVariable<Real>(tid, moose_obj->getParam<NonlinearVariableName>(NS::pressure))),
73 1466 : _enclosure_lm_var(moose_obj->isParamValid("enclosure_lm")
74 2060 : ? &sys.getScalarVariable(
75 2654 : tid, moose_obj->getParam<NonlinearVariableName>("enclosure_lm"))
76 : : nullptr),
77 : // dof indices
78 1466 : _qv_dof_indices(_grad_v_var.dofIndices()),
79 1466 : _v_dof_indices(_v_var.dofIndices()),
80 1466 : _lm_v_dof_indices(_v_face_var.dofIndices()),
81 1466 : _qw_dof_indices(_grad_w_var ? &_grad_w_var->dofIndices() : nullptr),
82 1466 : _w_dof_indices(_w_var ? &_w_var->dofIndices() : nullptr),
83 1466 : _lm_w_dof_indices(_w_face_var ? &_w_face_var->dofIndices() : nullptr),
84 1466 : _p_dof_indices(_pressure_var.dofIndices()),
85 1466 : _global_lm_dof_indices(_enclosure_lm_var ? &_enclosure_lm_var->dofIndices() : nullptr),
86 : // solutions
87 1466 : _qv_sol(_grad_v_var.sln()),
88 1466 : _v_sol(_v_var.sln()),
89 1466 : _lm_v_sol(_v_face_var.sln()),
90 1466 : _qw_sol(_grad_w_var ? &_grad_w_var->sln() : nullptr),
91 1466 : _w_sol(_w_var ? &_w_var->sln() : nullptr),
92 1466 : _lm_w_sol(_w_face_var ? &_w_face_var->sln() : nullptr),
93 1466 : _p_sol(_pressure_var.sln()),
94 2060 : _global_lm_dof_value(_enclosure_lm_var ? &_enclosure_lm_var->sln() : nullptr),
95 2932 : _rho(moose_obj->getParam<Real>(NS::density))
96 : {
97 1466 : if (mesh.dimension() > 2)
98 0 : mooseError("3D not yet implemented");
99 :
100 1466 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_v_var));
101 1466 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<RealVectorValue> &>(_grad_v_var));
102 1466 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_v_face_var));
103 1466 : mvdi->addMooseVariableDependency(&const_cast<MooseVariableFE<Real> &>(_pressure_var));
104 :
105 1466 : const auto vel_type = _u_var.feType();
106 2932 : auto check_type = [&vel_type](const auto & var)
107 : {
108 2932 : if (vel_type != var.feType())
109 0 : mooseError(
110 0 : 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 2932 : };
114 1466 : check_type(_v_var);
115 1466 : if (_w_var)
116 0 : check_type(*_w_var);
117 1466 : check_type(_pressure_var);
118 :
119 1466 : const auto grad_type = _grad_u_var.feType();
120 1466 : auto check_grad_type = [&grad_type](const auto & var)
121 : {
122 1466 : if (grad_type != var.feType())
123 0 : mooseError(
124 0 : 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 1466 : };
129 1466 : check_grad_type(_grad_v_var);
130 1466 : if (_grad_w_var)
131 0 : check_grad_type(*_grad_w_var);
132 :
133 1466 : const auto lm_type = _u_face_var.feType();
134 1466 : auto check_lm_type = [&lm_type](const auto & var)
135 : {
136 1466 : 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 1466 : };
142 1466 : check_lm_type(_v_face_var);
143 1466 : if (_w_face_var)
144 0 : check_lm_type(*_w_face_var);
145 1466 : }
146 :
147 : void
148 5786016 : 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 5786016 : DiffusionLHDGAssemblyHelper::scalarVolumeResidual(
158 : vel_gradient, body_force, JxW, qrule, current_elem, q_point, scalar_re);
159 :
160 28930080 : for (const auto qp : make_range(qrule.n_points()))
161 : {
162 23144064 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_u_sol, _v_sol, qp, vel_component);
163 : Gradient qp_p;
164 23144064 : qp_p(vel_component) = _p_sol[qp];
165 :
166 92576256 : for (const auto i : index_range(scalar_re))
167 : {
168 : // Scalar equation dependence on pressure dofs
169 69432192 : scalar_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * qp_p);
170 :
171 : // Scalar equation dependence on scalar dofs
172 69432192 : scalar_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * rho_vel_cross_vel);
173 : }
174 : }
175 5786016 : }
176 :
177 : void
178 1651840 : 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 1651840 : DiffusionLHDGAssemblyHelper::scalarVolumeJacobian(JxW, qrule, scalar_vector_jac);
187 :
188 6607360 : for (const auto i : make_range(scalar_vector_jac.m()))
189 24777600 : for (const auto qp : make_range(qrule.n_points()))
190 : {
191 : // Scalar equation dependence on pressure dofs
192 79288320 : for (const auto j : make_range(scalar_p_jac.n()))
193 : {
194 : Gradient p_phi;
195 59466240 : p_phi(vel_component) = _scalar_phi[j][qp];
196 59466240 : 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 79288320 : 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 59466240 : rhoVelCrossVelJacobian(_u_sol, _v_sol, qp, vel_component, 0, _scalar_phi, j);
207 59466240 : 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 59466240 : rhoVelCrossVelJacobian(_u_sol, _v_sol, qp, vel_component, 1, _scalar_phi, j);
213 59466240 : scalar_v_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * rho_vel_cross_vel);
214 : }
215 : }
216 : }
217 1651840 : }
218 :
219 : void
220 2893008 : 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 14465040 : for (const auto qp : make_range(qrule.n_points()))
230 : {
231 : // Prepare forcing function
232 11572032 : const auto f = pressure_mms_forcing_function(
233 11572032 : Moose::ElemQpArg{current_elem, qp, &qrule, q_point[qp]}, _ti.determineState());
234 :
235 11572032 : const Gradient vel(_u_sol[qp], _v_sol[qp]);
236 46288128 : for (const auto i : make_range(pressure_re.size()))
237 : {
238 34716096 : pressure_re(i) -= JxW[qp] * (_grad_scalar_phi[i][qp] * vel);
239 :
240 : // Pressure equation forcing function RHS
241 34716096 : pressure_re(i) -= JxW[qp] * _scalar_phi[i][qp] * f;
242 :
243 34716096 : 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 16677696 : pressure_re(i) -= JxW[qp] * _scalar_phi[i][qp] * (*_global_lm_dof_value)[0];
249 : }
250 : }
251 :
252 11572032 : if (_enclosure_lm_var)
253 5559232 : global_lm_re(0) -= JxW[qp] * _p_sol[qp];
254 : }
255 2893008 : }
256 :
257 : void
258 825920 : 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 4129600 : for (const auto qp : make_range(qrule.n_points()))
266 : {
267 13214720 : for (const auto i : make_range(p_u_vel_jac.m()))
268 : {
269 39644160 : for (const auto j : make_range(p_u_vel_jac.n()))
270 : {
271 : {
272 29733120 : const Gradient phi(_scalar_phi[j][qp], 0);
273 29733120 : 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 29733120 : p_v_vel_jac(i, j) -= JxW[qp] * (_grad_scalar_phi[i][qp] * phi);
278 : }
279 : }
280 9911040 : if (_enclosure_lm_var)
281 : p_global_lm_jac(i, 0) -= JxW[qp] * _scalar_phi[i][qp];
282 : }
283 :
284 3303680 : if (_enclosure_lm_var)
285 : {
286 5841920 : for (const auto j : make_range(global_lm_p_jac.n()))
287 4381440 : global_lm_p_jac(0, j) -= JxW[qp] * _scalar_phi[j][qp];
288 : }
289 : }
290 825920 : }
291 :
292 : RealVectorValue
293 91492672 : 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 91492672 : return _rho * U * U(vel_component);
300 : }
301 :
302 : RealVectorValue
303 1170955008 : 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 1170955008 : vector_phi(vel_j_component) = phi[j][qp];
314 : auto ret = vector_phi * U(vel_component);
315 1170955008 : if (vel_component == vel_j_component)
316 : ret += U * phi[j][qp];
317 : ret *= _rho;
318 1170955008 : return ret;
319 : }
320 :
321 : void
322 8543576 : NavierStokesLHDGAssemblyHelper::pressureFaceResidual(const MooseArray<Real> & JxW_face,
323 : const QBase & qrule_face,
324 : const MooseArray<Point> & normals,
325 : DenseVector<Number> & pressure_re)
326 : {
327 25630728 : for (const auto qp : make_range(qrule_face.n_points()))
328 : {
329 17087152 : const Gradient vel(_lm_u_sol[qp], _lm_v_sol[qp]);
330 : const auto vdotn = vel * normals[qp];
331 68348608 : for (const auto i : make_range(pressure_re.size()))
332 51261456 : pressure_re(i) += JxW_face[qp] * vdotn * _scalar_phi_face[i][qp];
333 : }
334 8543576 : }
335 :
336 : void
337 2438832 : 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 9755328 : for (const auto i : make_range(p_lm_u_vel_jac.m()))
348 51215472 : for (const auto j : make_range(p_lm_u_vel_jac.n()))
349 131696928 : for (const auto qp : make_range(qrule_face.n_points()))
350 : {
351 : {
352 87797952 : const Gradient phi(_lm_phi_face[j][qp], 0);
353 87797952 : 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 87797952 : p_lm_v_vel_jac(i, j) += JxW_face[qp] * phi * normals[qp] * _scalar_phi_face[i][qp];
358 : }
359 : }
360 2438832 : }
361 :
362 : void
363 17087152 : 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 17087152 : DiffusionLHDGAssemblyHelper::scalarFaceResidual(
373 : vector_sol, scalar_sol, lm_sol, JxW_face, qrule_face, normals, scalar_re);
374 :
375 51261456 : for (const auto qp : make_range(qrule_face.n_points()))
376 : {
377 : Gradient qp_p;
378 34174304 : qp_p(vel_component) = _p_sol[qp];
379 34174304 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_lm_u_sol, _lm_v_sol, qp, vel_component);
380 :
381 136697216 : for (const auto i : make_range(scalar_re.size()))
382 : {
383 : // pressure
384 102522912 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * (qp_p * normals[qp]);
385 :
386 : // lm from convection term
387 102522912 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
388 : }
389 : }
390 17087152 : }
391 :
392 : void
393 4877664 : 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 4877664 : DiffusionLHDGAssemblyHelper::scalarFaceJacobian(
406 : JxW_face, qrule_face, normals, scalar_vector_jac, scalar_scalar_jac, scalar_lm_jac);
407 :
408 19510656 : for (const auto i : make_range(scalar_lm_u_vel_jac.m()))
409 43898976 : for (const auto qp : make_range(qrule_face.n_points()))
410 : {
411 117063936 : for (const auto j : make_range(scalar_p_jac.n()))
412 : {
413 : Gradient p_phi;
414 87797952 : p_phi(vel_component) = _scalar_phi_face[j][qp];
415 : // pressure
416 87797952 : scalar_p_jac(i, j) += JxW_face[qp] * _scalar_phi_face[i][qp] * (p_phi * normals[qp]);
417 : }
418 :
419 204861888 : 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 175595904 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 0, _lm_phi_face, j);
429 175595904 : scalar_lm_u_vel_jac(i, j) +=
430 175595904 : 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 175595904 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 1, _lm_phi_face, j);
436 175595904 : scalar_lm_v_vel_jac(i, j) +=
437 175595904 : JxW_face[qp] * _scalar_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
438 : }
439 : }
440 : }
441 4877664 : }
442 :
443 : void
444 17087152 : 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 17087152 : DiffusionLHDGAssemblyHelper::lmFaceResidual(
455 : vector_sol, scalar_sol, lm_sol, JxW_face, qrule_face, normals, lm_re);
456 :
457 51261456 : for (const auto qp : make_range(qrule_face.n_points()))
458 : {
459 : Gradient qp_p;
460 34174304 : qp_p(vel_component) = _p_sol[qp];
461 34174304 : const auto rho_vel_cross_vel = rhoVelCrossVelResidual(_lm_u_sol, _lm_v_sol, qp, vel_component);
462 :
463 239220128 : for (const auto i : make_range(lm_re.size()))
464 : {
465 : // pressure
466 205045824 : 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 205045824 : if (neigh)
471 : // lm from convection term
472 204623232 : lm_re(i) += JxW_face[qp] * _lm_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
473 : }
474 : }
475 17087152 : }
476 :
477 : void
478 4877664 : 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 4877664 : DiffusionLHDGAssemblyHelper::lmFaceJacobian(
491 : JxW_face, qrule_face, normals, lm_vec_jac, lm_scalar_jac, lm_lm_jac);
492 :
493 34143648 : for (const auto i : make_range(lm_p_jac.m()))
494 87797952 : for (const auto qp : make_range(qrule_face.n_points()))
495 : {
496 234127872 : for (const auto j : make_range(lm_p_jac.n()))
497 : {
498 : Gradient p_phi;
499 175595904 : p_phi(vel_component) = _scalar_phi_face[j][qp];
500 175595904 : lm_p_jac(i, j) += JxW_face[qp] * _lm_phi_face[i][qp] * (p_phi * normals[qp]);
501 : }
502 :
503 409723776 : for (const auto j : make_range(lm_lm_u_vel_jac.n()))
504 351191808 : if (neigh)
505 : {
506 : // derivatives wrt 0th component of velocity
507 : {
508 : const auto rho_vel_cross_vel =
509 350415360 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 0, _lm_phi_face, j);
510 350415360 : lm_lm_u_vel_jac(i, j) +=
511 350415360 : 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 350415360 : rhoVelCrossVelJacobian(_lm_u_sol, _lm_v_sol, qp, vel_component, 1, _lm_phi_face, j);
517 350415360 : lm_lm_v_vel_jac(i, j) +=
518 350415360 : JxW_face[qp] * _lm_phi_face[i][qp] * rho_vel_cross_vel * normals[qp];
519 : }
520 : }
521 : }
522 4877664 : }
523 :
524 : void
525 135448 : 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 406344 : for (const auto qp : make_range(qrule_face.n_points()))
536 : {
537 : const Moose::ElemSideQpArg elem_side_qp_arg{
538 270896 : current_elem, current_side, qp, &qrule_face, q_point_face[qp]};
539 270896 : const auto time_arg = _ti.determineState();
540 270896 : const RealVectorValue dirichlet_velocity((*dirichlet_vel[0])(elem_side_qp_arg, time_arg),
541 270896 : (*dirichlet_vel[1])(elem_side_qp_arg, time_arg),
542 270896 : (*dirichlet_vel[2])(elem_side_qp_arg, time_arg));
543 : const auto vdotn = dirichlet_velocity * normals[qp];
544 1083584 : for (const auto i : make_range(pressure_re.size()))
545 812688 : pressure_re(i) += JxW_face[qp] * vdotn * _scalar_phi_face[i][qp];
546 : }
547 135448 : }
548 :
549 : void
550 270896 : 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 270896 : DiffusionLHDGAssemblyHelper::scalarDirichletResidual(vector_sol,
564 : scalar_sol,
565 270896 : *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 812688 : for (const auto qp : make_range(qrule_face.n_points()))
575 : {
576 : Gradient qp_p;
577 541792 : qp_p(vel_component) = _p_sol[qp];
578 :
579 : const Moose::ElemSideQpArg elem_side_qp_arg{
580 541792 : current_elem, current_side, qp, &qrule_face, q_point_face[qp]};
581 541792 : const auto time_arg = _ti.determineState();
582 541792 : const RealVectorValue dirichlet_velocity((*dirichlet_vel[0])(elem_side_qp_arg, time_arg),
583 541792 : (*dirichlet_vel[1])(elem_side_qp_arg, time_arg),
584 541792 : (*dirichlet_vel[2])(elem_side_qp_arg, time_arg));
585 541792 : const auto scalar_value = dirichlet_velocity(vel_component);
586 :
587 2167168 : for (const auto i : make_range(scalar_re.size()))
588 : {
589 : // pressure
590 1625376 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] * (qp_p * normals[qp]);
591 :
592 : // dirichlet lm from advection term
593 1625376 : scalar_re(i) += JxW_face[qp] * _scalar_phi_face[i][qp] *
594 1625376 : (_rho * dirichlet_velocity * normals[qp]) * scalar_value;
595 : }
596 : }
597 270896 : }
598 :
599 : void
600 77856 : 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 77856 : DiffusionLHDGAssemblyHelper::scalarDirichletJacobian(
609 : JxW_face, qrule_face, normals, scalar_vector_jac, scalar_scalar_jac);
610 :
611 311424 : for (const auto i : make_range(scalar_pressure_jac.m()))
612 934272 : for (const auto j : make_range(scalar_pressure_jac.n()))
613 2102112 : for (const auto qp : make_range(qrule_face.n_points()))
614 : {
615 : Gradient p_phi;
616 1401408 : p_phi(vel_component) = _scalar_phi_face[j][qp];
617 : // pressure
618 1401408 : scalar_pressure_jac(i, j) += JxW_face[qp] * _scalar_phi_face[i][qp] * (p_phi * normals[qp]);
619 : }
620 77856 : }
|
