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 "NumericalFlux3EqnDGKernel.h"
11 : #include "NumericalFlux3EqnBase.h"
12 : #include "MooseVariable.h"
13 : #include "THMIndicesVACE.h"
14 :
15 : registerMooseObject("ThermalHydraulicsApp", NumericalFlux3EqnDGKernel);
16 :
17 : InputParameters
18 0 : NumericalFlux3EqnDGKernel::validParams()
19 : {
20 0 : InputParameters params = DGKernel::validParams();
21 :
22 0 : params.addClassDescription(
23 : "Adds side fluxes for the 1-D, 1-phase, variable-area Euler equations");
24 :
25 0 : params.addRequiredCoupledVar("A_linear", "Cross-sectional area, linear");
26 0 : params.addRequiredCoupledVar("rhoA", "Conserved variable: rho*A");
27 0 : params.addRequiredCoupledVar("rhouA", "Conserved variable: rho*u*A");
28 0 : params.addRequiredCoupledVar("rhoEA", "Conserved variable: rho*E*A");
29 :
30 0 : params.addRequiredParam<UserObjectName>("numerical_flux", "Name of numerical flux user object");
31 :
32 0 : return params;
33 0 : }
34 :
35 0 : NumericalFlux3EqnDGKernel::NumericalFlux3EqnDGKernel(const InputParameters & parameters)
36 : : DGKernel(parameters),
37 :
38 0 : _A_elem(coupledValue("A_linear")),
39 0 : _A_neig(coupledNeighborValue("A_linear")),
40 0 : _rhoA1(getMaterialProperty<Real>("rhoA")),
41 0 : _rhouA1(getMaterialProperty<Real>("rhouA")),
42 0 : _rhoEA1(getMaterialProperty<Real>("rhoEA")),
43 0 : _p1(getMaterialProperty<Real>("p")),
44 0 : _rhoA2(getNeighborMaterialProperty<Real>("rhoA")),
45 0 : _rhouA2(getNeighborMaterialProperty<Real>("rhouA")),
46 0 : _rhoEA2(getNeighborMaterialProperty<Real>("rhoEA")),
47 0 : _p2(getNeighborMaterialProperty<Real>("p")),
48 0 : _numerical_flux(getUserObject<NumericalFlux3EqnBase>("numerical_flux")),
49 0 : _rhoA_var(coupled("rhoA")),
50 0 : _rhouA_var(coupled("rhouA")),
51 0 : _rhoEA_var(coupled("rhoEA")),
52 0 : _jmap(getIndexMapping()),
53 0 : _equation_index(_jmap.at(_var.number()))
54 : {
55 0 : }
56 :
57 : Real
58 0 : NumericalFlux3EqnDGKernel::computeQpResidual(Moose::DGResidualType type)
59 : {
60 : // construct the left and right solution vectors from the reconstructed solution
61 0 : std::vector<Real> U1 = {_rhoA1[_qp], _rhouA1[_qp], _rhoEA1[_qp], _A_elem[_qp]};
62 0 : std::vector<Real> U2 = {_rhoA2[_qp], _rhouA2[_qp], _rhoEA2[_qp], _A_neig[_qp]};
63 :
64 0 : const Real nLR_dot_d = _current_side * 2 - 1.0;
65 :
66 : const std::vector<Real> & flux_elem =
67 0 : _numerical_flux.getFlux(_current_side, _current_elem->id(), true, U1, U2, nLR_dot_d);
68 : const std::vector<Real> & flux_neig =
69 0 : _numerical_flux.getFlux(_current_side, _current_elem->id(), false, U1, U2, nLR_dot_d);
70 :
71 : Real re = 0.0;
72 0 : switch (type)
73 : {
74 0 : case Moose::Element:
75 0 : re = flux_elem[_equation_index] * _test[_i][_qp];
76 0 : break;
77 0 : case Moose::Neighbor:
78 0 : re = -flux_neig[_equation_index] * _test_neighbor[_i][_qp];
79 0 : break;
80 : }
81 0 : return re;
82 : }
83 :
84 : Real
85 0 : NumericalFlux3EqnDGKernel::computeQpJacobian(Moose::DGJacobianType type)
86 : {
87 0 : return computeQpOffDiagJacobian(type, _var.number());
88 : }
89 :
90 : Real
91 0 : NumericalFlux3EqnDGKernel::computeQpOffDiagJacobian(Moose::DGJacobianType type, unsigned int jvar)
92 : {
93 : // construct the left and right solution vectors from the cell-average solution
94 0 : std::vector<Real> U1 = {_rhoA1[_qp], _rhouA1[_qp], _rhoEA1[_qp], _A_elem[_qp]};
95 0 : std::vector<Real> U2 = {_rhoA2[_qp], _rhouA2[_qp], _rhoEA2[_qp], _A_neig[_qp]};
96 :
97 : // Temporary hack to allow libMesh fix. After the fix, this should be changed to:
98 : // const Real nLR_dot_d = _normals[_qp] * _direction[_qp];
99 0 : const Real nLR_dot_d = _current_side * 2 - 1.0;
100 :
101 0 : const DenseMatrix<Real> & dF1_dU1 = _numerical_flux.getJacobian(
102 0 : true, true, _current_side, _current_elem->id(), U1, U2, nLR_dot_d);
103 0 : const DenseMatrix<Real> & dF1_dU2 = _numerical_flux.getJacobian(
104 0 : true, false, _current_side, _current_elem->id(), U1, U2, nLR_dot_d);
105 0 : const DenseMatrix<Real> & dF2_dU1 = _numerical_flux.getJacobian(
106 0 : false, true, _current_side, _current_elem->id(), U1, U2, nLR_dot_d);
107 0 : const DenseMatrix<Real> & dF2_dU2 = _numerical_flux.getJacobian(
108 0 : false, false, _current_side, _current_elem->id(), U1, U2, nLR_dot_d);
109 :
110 : Real re = 0.0;
111 0 : switch (type)
112 : {
113 0 : case Moose::ElementElement:
114 0 : re = dF1_dU1(_equation_index, _jmap.at(jvar)) * _phi[_j][_qp] * _test[_i][_qp];
115 0 : break;
116 0 : case Moose::ElementNeighbor:
117 0 : re = dF1_dU2(_equation_index, _jmap.at(jvar)) * _phi_neighbor[_j][_qp] * _test[_i][_qp];
118 0 : break;
119 0 : case Moose::NeighborElement:
120 0 : re = -dF2_dU1(_equation_index, _jmap.at(jvar)) * _phi[_j][_qp] * _test_neighbor[_i][_qp];
121 0 : break;
122 0 : case Moose::NeighborNeighbor:
123 0 : re = -dF2_dU2(_equation_index, _jmap.at(jvar)) * _phi_neighbor[_j][_qp] *
124 0 : _test_neighbor[_i][_qp];
125 0 : break;
126 : }
127 0 : return re;
128 : }
129 :
130 : std::map<unsigned int, unsigned int>
131 0 : NumericalFlux3EqnDGKernel::getIndexMapping() const
132 : {
133 : std::map<unsigned int, unsigned int> jmap;
134 0 : jmap.insert(std::pair<unsigned int, unsigned int>(_rhoA_var, THMVACE1D::MASS));
135 0 : jmap.insert(std::pair<unsigned int, unsigned int>(_rhouA_var, THMVACE1D::MOMENTUM));
136 0 : jmap.insert(std::pair<unsigned int, unsigned int>(_rhoEA_var, THMVACE1D::ENERGY));
137 :
138 0 : return jmap;
139 : }
|
