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 "ADGateValve1PhaseUserObject.h"
11 : #include "THMIndicesVACE.h"
12 : #include "ADNumericalFlux3EqnBase.h"
13 : #include "Function.h"
14 : #include "Numerics.h"
15 : #include "metaphysicl/parallel_numberarray.h"
16 : #include "metaphysicl/parallel_dualnumber.h"
17 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
18 : #include "libmesh/parallel_algebra.h"
19 :
20 : registerMooseObject("ThermalHydraulicsApp", ADGateValve1PhaseUserObject);
21 :
22 : const std::vector<std::pair<std::string, unsigned int>>
23 : ADGateValve1PhaseUserObject::_varname_eq_index_pairs{
24 : std::pair<std::string, unsigned int>("rhoA", THMVACE1D::MASS),
25 : std::pair<std::string, unsigned int>("rhouA", THMVACE1D::MOMENTUM),
26 : std::pair<std::string, unsigned int>("rhoEA", THMVACE1D::ENERGY)};
27 :
28 : InputParameters
29 166 : ADGateValve1PhaseUserObject::validParams()
30 : {
31 166 : InputParameters params = ADFlowJunctionUserObject::validParams();
32 :
33 332 : params.addRequiredParam<Real>("open_area_fraction",
34 : "Fraction of possible flow area that is open");
35 332 : params.addParam<Real>("open_area_fraction_min", 1e-10, "Minimum open area fraction");
36 :
37 332 : params.addRequiredCoupledVar("A", "Cross-sectional area of connected flow channels");
38 332 : params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
39 332 : params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
40 332 : params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
41 :
42 332 : params.addRequiredParam<UserObjectName>("numerical_flux", "Numerical flux user object name");
43 :
44 332 : params.addRequiredParam<std::string>("component_name", "Name of the associated component");
45 :
46 332 : params.addClassDescription("Gate valve user object for 1-phase flow");
47 :
48 332 : params.declareControllable("open_area_fraction");
49 :
50 166 : return params;
51 0 : }
52 :
53 90 : ADGateValve1PhaseUserObject::ADGateValve1PhaseUserObject(const InputParameters & params)
54 : : ADFlowJunctionUserObject(params),
55 :
56 90 : _f_open(getParam<Real>("open_area_fraction")),
57 180 : _f_open_min(getParam<Real>("open_area_fraction_min")),
58 :
59 90 : _A(adCoupledValue("A")),
60 90 : _rhoA(adCoupledValue("rhoA")),
61 90 : _rhouA(adCoupledValue("rhouA")),
62 90 : _rhoEA(adCoupledValue("rhoEA")),
63 :
64 90 : _rhoA_jvar(coupled("rhoA")),
65 90 : _rhouA_jvar(coupled("rhouA")),
66 90 : _rhoEA_jvar(coupled("rhoEA")),
67 :
68 180 : _p(getADMaterialProperty<Real>("p")),
69 :
70 90 : _numerical_flux(getUserObject<ADNumericalFlux3EqnBase>("numerical_flux")),
71 :
72 180 : _component_name(getParam<std::string>("component_name")),
73 :
74 180 : _dir(getMaterialProperty<RealVectorValue>("direction")),
75 :
76 90 : _solutions(_n_connections),
77 90 : _fluxes(_n_connections, std::vector<ADReal>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
78 90 : _dof_indices(_n_connections, std::vector<dof_id_type>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
79 :
80 90 : _stored_p(_n_connections),
81 :
82 90 : _elem_ids(_n_connections),
83 90 : _local_side_ids(_n_connections),
84 :
85 90 : _areas(_n_connections),
86 180 : _directions(_n_connections)
87 : {
88 90 : }
89 :
90 : void
91 10382 : ADGateValve1PhaseUserObject::initialize()
92 : {
93 : _connection_indices.clear();
94 :
95 10382 : std::vector<ADReal> zero(THMVACE1D::N_FLUX_INPUTS, ADReal(0.));
96 31146 : for (auto & s : _solutions)
97 20764 : s = zero;
98 41528 : for (unsigned int i = 0; i < THMVACE1D::N_FLUX_OUTPUTS; i++)
99 : {
100 31146 : _fluxes[0][i] = 0;
101 31146 : _fluxes[1][i] = 0;
102 : }
103 10382 : }
104 :
105 : void
106 9892 : ADGateValve1PhaseUserObject::execute()
107 : {
108 : // Get connection index
109 9892 : const unsigned int c = getBoundaryIDIndex();
110 9892 : _connection_indices.push_back(c);
111 :
112 : // Store DoF indices
113 39568 : for (const auto & varname_eq_index_pair : _varname_eq_index_pairs)
114 : {
115 29676 : MooseVariable * var = getVar(varname_eq_index_pair.first, 0);
116 : auto && dofs = var->dofIndices();
117 : mooseAssert(dofs.size() == 1, "There should be only one DoF index.");
118 29676 : _dof_indices[c][varname_eq_index_pair.second] = dofs[0];
119 : }
120 :
121 : // Store solution vector for connection
122 9892 : std::vector<ADReal> U(THMVACE1D::N_FLUX_INPUTS, ADReal(0.));
123 9892 : U[THMVACE1D::RHOA] = _rhoA[0];
124 9892 : U[THMVACE1D::RHOUA] = _rhouA[0];
125 9892 : U[THMVACE1D::RHOEA] = _rhoEA[0];
126 9892 : U[THMVACE1D::AREA] = _A[0];
127 9892 : _solutions[c] = U;
128 :
129 : // Store element ID and local side ID for connection
130 9892 : _elem_ids[c] = _current_elem->id();
131 9892 : _local_side_ids[c] = _current_side;
132 :
133 : // Store direction and area of channel at junction (used for error-checking)
134 9892 : _areas[c] = _A[0];
135 9892 : _directions[c] = _dir[0];
136 :
137 9892 : _stored_p[c] = _p[0];
138 9892 : }
139 :
140 : void
141 1934 : ADGateValve1PhaseUserObject::threadJoin(const UserObject & uo)
142 : {
143 : const auto & junction_uo = static_cast<const ADGateValve1PhaseUserObject &>(uo);
144 :
145 : // Store the data computed/retrieved in the other threads
146 4109 : for (unsigned int i = 0; i < junction_uo._connection_indices.size(); i++)
147 : {
148 2175 : const unsigned int c = junction_uo._connection_indices[i];
149 :
150 2175 : _dof_indices[c] = junction_uo._dof_indices[c];
151 2175 : _solutions[c] = junction_uo._solutions[c];
152 2175 : _elem_ids[c] = junction_uo._elem_ids[c];
153 2175 : _local_side_ids[c] = junction_uo._local_side_ids[c];
154 2175 : _areas[c] = junction_uo._areas[c];
155 2175 : _directions[c] = junction_uo._directions[c];
156 2175 : _stored_p[c] = junction_uo._stored_p[c];
157 : }
158 1934 : }
159 :
160 : void
161 8448 : ADGateValve1PhaseUserObject::finalize()
162 : {
163 : // Check direction compatibility
164 8448 : if (!THM::areParallelVectors(_directions[0], _directions[1]))
165 2 : mooseError(_component_name, ": The connected channels must be parallel at the junction.");
166 :
167 25338 : for (unsigned int i = 0; i < _n_connections; i++)
168 : {
169 16892 : processor_id_type owner_proc = _processor_ids[i];
170 16892 : comm().broadcast(_elem_ids[i], owner_proc, true);
171 16892 : comm().broadcast(_local_side_ids[i], owner_proc, true);
172 16892 : comm().broadcast(_solutions[i], owner_proc, true);
173 16892 : comm().broadcast(_directions[i], owner_proc, true);
174 16892 : comm().broadcast(_areas[i], owner_proc, true);
175 16892 : comm().broadcast(_stored_p[i], owner_proc, true);
176 : }
177 :
178 8446 : const Real & n1_dot_d1 = _normal[0];
179 8446 : const Real d1_dot_d2 = _directions[0] * _directions[1];
180 :
181 : const ADReal & A1 = _areas[0];
182 : const ADReal & A2 = _areas[1];
183 8446 : const ADReal A = std::min(A1, A2);
184 8446 : const ADReal A_flow = _f_open * A;
185 :
186 8446 : if (_f_open > _f_open_min)
187 : {
188 : // compute flow contribution
189 4196 : std::vector<ADReal> U_flow1 = _solutions[0];
190 4196 : std::vector<ADReal> U_flow2 = _solutions[1];
191 4196 : U_flow1[THMVACE1D::AREA] = A_flow;
192 4196 : U_flow2[THMVACE1D::AREA] = A_flow;
193 16784 : for (unsigned int i = 0; i < THMVACE1D::N_FLUX_OUTPUTS; i++)
194 : {
195 12588 : U_flow1[i] *= A_flow / A1;
196 12588 : U_flow2[i] *= A_flow / A2;
197 : }
198 4196 : U_flow2[THMVACE1D::RHOUA] *= d1_dot_d2;
199 :
200 4196 : _fluxes[0] = _numerical_flux.getFlux(
201 4196 : _local_side_ids[0], _elem_ids[0], true, U_flow1, U_flow2, n1_dot_d1);
202 :
203 4196 : _fluxes[1] = _numerical_flux.getFlux(
204 4196 : _local_side_ids[0], _elem_ids[0], false, U_flow1, U_flow2, n1_dot_d1);
205 4196 : _fluxes[1][THMVACE1D::RHOA] *= d1_dot_d2;
206 4196 : _fluxes[1][THMVACE1D::RHOEA] *= d1_dot_d2;
207 : }
208 :
209 : // compute wall contribution
210 25338 : for (unsigned int c = 0; c < _n_connections; c++)
211 : {
212 16892 : const ADReal A_wall = _areas[c] - A_flow;
213 16892 : _fluxes[c][THMVACE1D::RHOUA] += _stored_p[c] * A_wall;
214 : }
215 8446 : }
216 :
217 : const std::vector<ADReal> &
218 29460 : ADGateValve1PhaseUserObject::getFlux(const unsigned int & connection_index) const
219 : {
220 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
221 :
222 29460 : return _fluxes[connection_index];
223 : }
|
