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 "ADJunctionOneToOne1PhaseUserObject.h"
11 : #include "THMIndicesVACE.h"
12 : #include "FlowModel1PhaseUtils.h"
13 : #include "SinglePhaseFluidProperties.h"
14 : #include "ADNumericalFlux3EqnBase.h"
15 : #include "Numerics.h"
16 : #include "metaphysicl/parallel_numberarray.h"
17 : #include "metaphysicl/parallel_dualnumber.h"
18 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
19 : #include "libmesh/parallel_algebra.h"
20 :
21 : registerMooseObject("ThermalHydraulicsApp", ADJunctionOneToOne1PhaseUserObject);
22 :
23 : const std::vector<std::pair<std::string, unsigned int>>
24 : ADJunctionOneToOne1PhaseUserObject::_varname_eq_index_pairs{
25 : std::pair<std::string, unsigned int>("rhoA", THMVACE1D::MASS),
26 : std::pair<std::string, unsigned int>("rhouA", THMVACE1D::MOMENTUM),
27 : std::pair<std::string, unsigned int>("rhoEA", THMVACE1D::ENERGY)};
28 :
29 : Threads::spin_mutex ADJunctionOneToOne1PhaseUserObject::_spin_mutex;
30 :
31 : InputParameters
32 219 : ADJunctionOneToOne1PhaseUserObject::validParams()
33 : {
34 219 : InputParameters params = ADFlowJunctionUserObject::validParams();
35 219 : params += SlopeReconstruction1DInterface<true>::validParams();
36 :
37 438 : params.addRequiredCoupledVar(
38 : "A_elem", "Piecewise constant cross-sectional area of connected flow channels");
39 438 : params.addRequiredCoupledVar("A_linear",
40 : "Piecewise linear cross-sectional area of connected flow channels");
41 438 : params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
42 438 : params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
43 438 : params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
44 :
45 438 : params.addRequiredParam<UserObjectName>("fluid_properties",
46 : "Name of fluid properties user object");
47 :
48 438 : params.addRequiredParam<UserObjectName>("numerical_flux", "Numerical flux user object name");
49 :
50 438 : params.addRequiredParam<std::string>("junction_name", "Name of the junction component");
51 :
52 219 : params.addClassDescription(
53 : "Computes flux between two subdomains for 1-phase one-to-one junction");
54 :
55 219 : return params;
56 0 : }
57 :
58 116 : ADJunctionOneToOne1PhaseUserObject::ADJunctionOneToOne1PhaseUserObject(
59 116 : const InputParameters & params)
60 : : ADFlowJunctionUserObject(params),
61 : SlopeReconstruction1DInterface<true>(this),
62 :
63 116 : _A_linear(adCoupledValue("A_linear")),
64 :
65 116 : _A_var(getVar("A_elem", 0)),
66 116 : _rhoA_var(getVar("rhoA", 0)),
67 116 : _rhouA_var(getVar("rhouA", 0)),
68 116 : _rhoEA_var(getVar("rhoEA", 0)),
69 :
70 116 : _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties")),
71 :
72 116 : _numerical_flux(getUserObject<ADNumericalFlux3EqnBase>("numerical_flux")),
73 :
74 232 : _junction_name(getParam<std::string>("junction_name")),
75 :
76 232 : _dir(getMaterialProperty<RealVectorValue>("direction")),
77 :
78 116 : _primitive_solutions(_n_connections),
79 116 : _neighbor_primitive_solutions(_n_connections),
80 :
81 116 : _fluxes(_n_connections),
82 116 : _dof_indices(_n_connections, std::vector<dof_id_type>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
83 :
84 116 : _elem_ids(_n_connections),
85 116 : _local_side_ids(_n_connections),
86 :
87 116 : _areas_linear(_n_connections),
88 116 : _directions(_n_connections),
89 116 : _positions(_n_connections),
90 116 : _neighbor_positions(_n_connections),
91 116 : _delta_x(_n_connections),
92 232 : _has_neighbor(_n_connections)
93 : {
94 116 : _U_vars.resize(THMVACE1D::N_FLUX_INPUTS);
95 116 : _U_vars[THMVACE1D::RHOA] = _rhoA_var;
96 116 : _U_vars[THMVACE1D::RHOUA] = _rhouA_var;
97 116 : _U_vars[THMVACE1D::RHOEA] = _rhoEA_var;
98 116 : _U_vars[THMVACE1D::AREA] = _A_var;
99 116 : }
100 :
101 : void
102 5670 : ADJunctionOneToOne1PhaseUserObject::initialize()
103 : {
104 5670 : _connection_indices.clear();
105 :
106 : // Broadcasts require vectors on all processes to have the same size
107 5670 : std::vector<ADReal> zero(THMVACE1D::N_PRIM_VARS, ADReal(0.));
108 17010 : for (unsigned int c = 0; c < _n_connections; c++)
109 : {
110 11340 : _primitive_solutions[c] = zero;
111 11340 : _neighbor_primitive_solutions[c] = zero;
112 : }
113 5670 : }
114 :
115 : void
116 8512 : ADJunctionOneToOne1PhaseUserObject::execute()
117 : {
118 : // Get connection index
119 8512 : const unsigned int c = getBoundaryIDIndex();
120 8512 : _connection_indices.push_back(c);
121 :
122 : // Store DoF indices
123 34048 : for (const auto & varname_eq_index_pair : _varname_eq_index_pairs)
124 : {
125 25536 : MooseVariable * var = getVar(varname_eq_index_pair.first, 0);
126 : auto && dofs = var->dofIndices();
127 : mooseAssert(dofs.size() == 1, "There should be only one DoF index.");
128 25536 : _dof_indices[c][varname_eq_index_pair.second] = dofs[0];
129 : }
130 :
131 : // Store primitive solution vectors for connection
132 : const auto U_avg =
133 8512 : FlowModel1PhaseUtils::getElementalSolutionVector<true>(_current_elem, _U_vars, _is_implicit);
134 8512 : _primitive_solutions[c] = FlowModel1PhaseUtils::computePrimitiveSolutionVector<true>(U_avg, _fp);
135 :
136 : // Get the neighbor solution and position. There should be one neighbor, unless
137 : // there is only one element in the subdomain. Because broadcasting requires
138 : // data sizes to be allocated on all processes, we cannot work with a vector
139 : // of a size not known upfront, so we work with a single neighbor and have
140 : // a flag that states whether there is one or zero neighbors.
141 : std::vector<std::vector<GenericReal<true>>> W_neighbor;
142 : std::vector<Point> x_neighbor;
143 8512 : getNeighborPrimitiveVariables(_current_elem, W_neighbor, x_neighbor);
144 8512 : if (W_neighbor.size() == 1)
145 : {
146 8480 : _neighbor_primitive_solutions[c] = W_neighbor[0];
147 8480 : _neighbor_positions[c] = x_neighbor[0];
148 : _has_neighbor[c] = true;
149 : }
150 : else
151 : _has_neighbor[c] = false;
152 :
153 : // Store element ID and local side ID for connection
154 8512 : _elem_ids[c] = _current_elem->id();
155 8512 : _local_side_ids[c] = _current_side;
156 :
157 : // Store direction, area, and position of channel at junction
158 8512 : _areas_linear[c] = _A_linear[0];
159 8512 : _directions[c] = _dir[0];
160 8512 : _positions[c] = _q_point[0];
161 8512 : _delta_x[c] = (_positions[c] - _current_elem->vertex_average()) * _directions[c];
162 8512 : }
163 :
164 : void
165 741 : ADJunctionOneToOne1PhaseUserObject::threadJoin(const UserObject & uo)
166 : {
167 : const auto & junction_uo = static_cast<const ADJunctionOneToOne1PhaseUserObject &>(uo);
168 :
169 : // Store the data computed/retrieved in the other threads
170 1480 : for (unsigned int i = 0; i < junction_uo._connection_indices.size(); i++)
171 : {
172 739 : const unsigned int c = junction_uo._connection_indices[i];
173 :
174 739 : _dof_indices[c] = junction_uo._dof_indices[c];
175 739 : _primitive_solutions[c] = junction_uo._primitive_solutions[c];
176 739 : _neighbor_primitive_solutions[c] = junction_uo._neighbor_primitive_solutions[c];
177 739 : _elem_ids[c] = junction_uo._elem_ids[c];
178 739 : _local_side_ids[c] = junction_uo._local_side_ids[c];
179 739 : _areas_linear[c] = junction_uo._areas_linear[c];
180 739 : _directions[c] = junction_uo._directions[c];
181 739 : _positions[c] = junction_uo._positions[c];
182 739 : _neighbor_positions[c] = junction_uo._neighbor_positions[c];
183 739 : _delta_x[c] = junction_uo._delta_x[c];
184 739 : _has_neighbor[c] = junction_uo._has_neighbor[c];
185 : }
186 741 : }
187 :
188 : void
189 4929 : ADJunctionOneToOne1PhaseUserObject::finalize()
190 : {
191 14787 : for (unsigned int i = 0; i < _n_connections; i++)
192 : {
193 9858 : processor_id_type owner_proc = _processor_ids[i];
194 9858 : comm().broadcast(_primitive_solutions[i], owner_proc, true);
195 9858 : comm().broadcast(_neighbor_primitive_solutions[i], owner_proc, true);
196 9858 : comm().broadcast(_elem_ids[i], owner_proc, true);
197 9858 : comm().broadcast(_local_side_ids[i], owner_proc, true);
198 9858 : comm().broadcast(_areas_linear[i], owner_proc, true);
199 9858 : comm().broadcast(_directions[i], owner_proc, true);
200 9858 : comm().broadcast(_positions[i], owner_proc, true);
201 9858 : comm().broadcast(_neighbor_positions[i], owner_proc, true);
202 9858 : comm().broadcast(_delta_x[i], owner_proc, true);
203 :
204 : // Vectors of bools are special, so we must broadcast a single bool instead
205 9858 : bool has_neighbor = _has_neighbor[i];
206 9858 : comm().broadcast(has_neighbor, owner_proc, true);
207 9858 : _has_neighbor[i] = has_neighbor;
208 : }
209 :
210 : const auto & W0_avg = _primitive_solutions[0];
211 : const auto & W1_avg = _primitive_solutions[1];
212 :
213 : // Multiplier for possibly different coordinate systems
214 4929 : const Real dir_mult = -_normal[0] * _normal[1];
215 4929 : auto W0_ref = W0_avg;
216 4929 : auto W1_ref = W1_avg;
217 4929 : W0_ref[THMVACE1D::VELOCITY] *= dir_mult;
218 4929 : W1_ref[THMVACE1D::VELOCITY] *= dir_mult;
219 :
220 : std::vector<std::vector<GenericReal<true>>> W_neighbor0;
221 : std::vector<Point> x_neighbor0;
222 4929 : if (_has_neighbor[0])
223 : {
224 4905 : W_neighbor0.push_back(_neighbor_primitive_solutions[0]);
225 4905 : x_neighbor0.push_back(_neighbor_positions[0]);
226 : }
227 :
228 : std::vector<std::vector<GenericReal<true>>> W_neighbor1;
229 : std::vector<Point> x_neighbor1;
230 4929 : if (_has_neighbor[1])
231 : {
232 4905 : W_neighbor1.push_back(_neighbor_primitive_solutions[1]);
233 4905 : x_neighbor1.push_back(_neighbor_positions[1]);
234 : }
235 :
236 : const auto slopes0 = getBoundaryElementSlopes(
237 4929 : W0_avg, _positions[0], _directions[0], W_neighbor0, x_neighbor0, W1_ref);
238 : const auto slopes1 = getBoundaryElementSlopes(
239 4929 : W1_avg, _positions[1], _directions[1], W_neighbor1, x_neighbor1, W0_ref);
240 :
241 4929 : auto W0 = W0_avg;
242 4929 : auto W1 = W1_avg;
243 19716 : for (unsigned int m = 0; m < THMVACE1D::N_PRIM_VARS; m++)
244 : {
245 29574 : W0[m] = W0_avg[m] + slopes0[m] * _delta_x[0];
246 14787 : W1[m] = W1_avg[m] + slopes1[m] * _delta_x[1];
247 : }
248 :
249 : auto U0 =
250 4929 : FlowModel1PhaseUtils::computeConservativeSolutionVector<true>(W0, _areas_linear[0], _fp);
251 : auto U1 =
252 4929 : FlowModel1PhaseUtils::computeConservativeSolutionVector<true>(W1, _areas_linear[1], _fp);
253 4929 : U1[THMVACE1D::RHOUA] *= dir_mult;
254 :
255 4929 : _fluxes[0] = _numerical_flux.getFlux(_local_side_ids[0], _elem_ids[0], true, U0, U1, _normal[0]);
256 :
257 4929 : _fluxes[1] = _numerical_flux.getFlux(_local_side_ids[0], _elem_ids[0], false, U0, U1, _normal[0]);
258 4929 : _fluxes[1][THMVACE1D::RHOA] *= dir_mult;
259 4929 : _fluxes[1][THMVACE1D::RHOEA] *= dir_mult;
260 4929 : }
261 :
262 : const std::vector<ADReal> &
263 25122 : ADJunctionOneToOne1PhaseUserObject::getFlux(const unsigned int & connection_index) const
264 : {
265 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
266 :
267 25122 : return _fluxes[connection_index];
268 : }
269 :
270 : std::vector<ADReal>
271 8480 : ADJunctionOneToOne1PhaseUserObject::computeElementPrimitiveVariables(const Elem * elem) const
272 : {
273 : const auto U =
274 8480 : FlowModel1PhaseUtils::getElementalSolutionVector<true>(elem, _U_vars, _is_implicit);
275 16960 : return FlowModel1PhaseUtils::computePrimitiveSolutionVector<true>(U, _fp);
276 8480 : }
|
