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 795 : ADJunctionOneToOne1PhaseUserObject::validParams()
33 : {
34 795 : InputParameters params = ADFlowJunctionUserObject::validParams();
35 795 : params += SlopeReconstruction1DInterface<true>::validParams();
36 :
37 1590 : params.addRequiredCoupledVar(
38 : "A_elem", "Piecewise constant cross-sectional area of connected flow channels");
39 1590 : params.addRequiredCoupledVar("A_linear",
40 : "Piecewise linear cross-sectional area of connected flow channels");
41 1590 : params.addRequiredCoupledVar("rhoA", "rho*A of the connected flow channels");
42 1590 : params.addRequiredCoupledVar("rhouA", "rhou*A of the connected flow channels");
43 1590 : params.addRequiredCoupledVar("rhoEA", "rhoE*A of the connected flow channels");
44 :
45 1590 : params.addRequiredParam<UserObjectName>("fluid_properties",
46 : "Name of fluid properties user object");
47 :
48 1590 : params.addRequiredParam<UserObjectName>("numerical_flux", "Numerical flux user object name");
49 :
50 1590 : params.addRequiredParam<std::string>("junction_name", "Name of the junction component");
51 :
52 795 : params.addClassDescription(
53 : "Computes flux between two subdomains for 1-phase one-to-one junction");
54 :
55 795 : return params;
56 0 : }
57 :
58 434 : ADJunctionOneToOne1PhaseUserObject::ADJunctionOneToOne1PhaseUserObject(
59 434 : const InputParameters & params)
60 : : ADFlowJunctionUserObject(params),
61 : SlopeReconstruction1DInterface<true>(this),
62 :
63 434 : _A_linear(adCoupledValue("A_linear")),
64 :
65 434 : _A_var(getVar("A_elem", 0)),
66 434 : _rhoA_var(getVar("rhoA", 0)),
67 434 : _rhouA_var(getVar("rhouA", 0)),
68 434 : _rhoEA_var(getVar("rhoEA", 0)),
69 :
70 434 : _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties")),
71 :
72 434 : _numerical_flux(getUserObject<ADNumericalFlux3EqnBase>("numerical_flux")),
73 :
74 868 : _junction_name(getParam<std::string>("junction_name")),
75 :
76 868 : _dir(getMaterialProperty<RealVectorValue>("direction")),
77 :
78 434 : _primitive_solutions(_n_connections),
79 434 : _neighbor_primitive_solutions(_n_connections),
80 :
81 434 : _fluxes(_n_connections),
82 434 : _dof_indices(_n_connections, std::vector<dof_id_type>(THMVACE1D::N_FLUX_OUTPUTS, 0)),
83 :
84 434 : _elem_ids(_n_connections),
85 434 : _local_side_ids(_n_connections),
86 :
87 434 : _areas_linear(_n_connections),
88 434 : _directions(_n_connections),
89 434 : _positions(_n_connections),
90 434 : _neighbor_positions(_n_connections),
91 434 : _delta_x(_n_connections),
92 1302 : _has_neighbor(_n_connections)
93 : {
94 434 : _U_vars.resize(THMVACE1D::N_FLUX_INPUTS);
95 434 : _U_vars[THMVACE1D::RHOA] = _rhoA_var;
96 434 : _U_vars[THMVACE1D::RHOUA] = _rhouA_var;
97 434 : _U_vars[THMVACE1D::RHOEA] = _rhoEA_var;
98 434 : _U_vars[THMVACE1D::AREA] = _A_var;
99 434 : }
100 :
101 : void
102 38772 : ADJunctionOneToOne1PhaseUserObject::initialize()
103 : {
104 : _connection_indices.clear();
105 :
106 : // Broadcasts require vectors on all processes to have the same size
107 38772 : std::vector<ADReal> zero(THMVACE1D::N_PRIM_VARS, ADReal(0.));
108 116316 : for (unsigned int c = 0; c < _n_connections; c++)
109 : {
110 77544 : _primitive_solutions[c] = zero;
111 77544 : _neighbor_primitive_solutions[c] = zero;
112 : }
113 38772 : }
114 :
115 : void
116 56126 : ADJunctionOneToOne1PhaseUserObject::execute()
117 : {
118 : // Get connection index
119 56126 : const unsigned int c = getBoundaryIDIndex();
120 56126 : _connection_indices.push_back(c);
121 :
122 : // Store DoF indices
123 224504 : for (const auto & varname_eq_index_pair : _varname_eq_index_pairs)
124 : {
125 168378 : 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 168378 : _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 56126 : FlowModel1PhaseUtils::getElementalSolutionVector<true>(_current_elem, _U_vars, _is_implicit);
134 112252 : _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 56126 : getNeighborPrimitiveVariables(_current_elem, W_neighbor, x_neighbor);
144 56126 : if (W_neighbor.size() == 1)
145 : {
146 56086 : _neighbor_primitive_solutions[c] = W_neighbor[0];
147 56086 : _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 56126 : _elem_ids[c] = _current_elem->id();
155 56126 : _local_side_ids[c] = _current_side;
156 :
157 : // Store direction, area, and position of channel at junction
158 56126 : _areas_linear[c] = _A_linear[0];
159 56126 : _directions[c] = _dir[0];
160 56126 : _positions[c] = _q_point[0];
161 56126 : _delta_x[c] = (_positions[c] - _current_elem->vertex_average()) * _directions[c];
162 112252 : }
163 :
164 : void
165 6221 : 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 12733 : for (unsigned int i = 0; i < junction_uo._connection_indices.size(); i++)
171 : {
172 6512 : const unsigned int c = junction_uo._connection_indices[i];
173 :
174 6512 : _dof_indices[c] = junction_uo._dof_indices[c];
175 6512 : _primitive_solutions[c] = junction_uo._primitive_solutions[c];
176 6512 : _neighbor_primitive_solutions[c] = junction_uo._neighbor_primitive_solutions[c];
177 6512 : _elem_ids[c] = junction_uo._elem_ids[c];
178 6512 : _local_side_ids[c] = junction_uo._local_side_ids[c];
179 6512 : _areas_linear[c] = junction_uo._areas_linear[c];
180 6512 : _directions[c] = junction_uo._directions[c];
181 6512 : _positions[c] = junction_uo._positions[c];
182 6512 : _neighbor_positions[c] = junction_uo._neighbor_positions[c];
183 6512 : _delta_x[c] = junction_uo._delta_x[c];
184 6512 : _has_neighbor[c] = junction_uo._has_neighbor[c];
185 : }
186 6221 : }
187 :
188 : void
189 32551 : ADJunctionOneToOne1PhaseUserObject::finalize()
190 : {
191 97653 : for (unsigned int i = 0; i < _n_connections; i++)
192 : {
193 65102 : processor_id_type owner_proc = _processor_ids[i];
194 65102 : comm().broadcast(_primitive_solutions[i], owner_proc, true);
195 65102 : comm().broadcast(_neighbor_primitive_solutions[i], owner_proc, true);
196 65102 : comm().broadcast(_elem_ids[i], owner_proc, true);
197 65102 : comm().broadcast(_local_side_ids[i], owner_proc, true);
198 65102 : comm().broadcast(_areas_linear[i], owner_proc, true);
199 65102 : comm().broadcast(_directions[i], owner_proc, true);
200 65102 : comm().broadcast(_positions[i], owner_proc, true);
201 65102 : comm().broadcast(_neighbor_positions[i], owner_proc, true);
202 65102 : comm().broadcast(_delta_x[i], owner_proc, true);
203 :
204 : // Vectors of bools are special, so we must broadcast a single bool instead
205 65102 : bool has_neighbor = _has_neighbor[i];
206 65102 : comm().broadcast(has_neighbor, owner_proc, true);
207 65102 : _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 32551 : const Real dir_mult = -_normal[0] * _normal[1];
215 32551 : auto W0_ref = W0_avg;
216 32551 : auto W1_ref = W1_avg;
217 32551 : W0_ref[THMVACE1D::VELOCITY] *= dir_mult;
218 32551 : W1_ref[THMVACE1D::VELOCITY] *= dir_mult;
219 :
220 : std::vector<std::vector<GenericReal<true>>> W_neighbor0;
221 : std::vector<Point> x_neighbor0;
222 32551 : if (_has_neighbor[0])
223 : {
224 32519 : W_neighbor0.push_back(_neighbor_primitive_solutions[0]);
225 32519 : 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 32551 : if (_has_neighbor[1])
231 : {
232 32519 : W_neighbor1.push_back(_neighbor_primitive_solutions[1]);
233 32519 : x_neighbor1.push_back(_neighbor_positions[1]);
234 : }
235 :
236 : const auto slopes0 = getBoundaryElementSlopes(
237 32551 : W0_avg, _positions[0], _directions[0], W_neighbor0, x_neighbor0, W1_ref);
238 : const auto slopes1 = getBoundaryElementSlopes(
239 32551 : W1_avg, _positions[1], _directions[1], W_neighbor1, x_neighbor1, W0_ref);
240 :
241 32551 : auto W0 = W0_avg;
242 32551 : auto W1 = W1_avg;
243 130204 : for (unsigned int m = 0; m < THMVACE1D::N_PRIM_VARS; m++)
244 : {
245 195306 : W0[m] = W0_avg[m] + slopes0[m] * _delta_x[0];
246 97653 : W1[m] = W1_avg[m] + slopes1[m] * _delta_x[1];
247 : }
248 :
249 : auto U0 =
250 32551 : FlowModel1PhaseUtils::computeConservativeSolutionVector<true>(W0, _areas_linear[0], _fp);
251 : auto U1 =
252 32551 : FlowModel1PhaseUtils::computeConservativeSolutionVector<true>(W1, _areas_linear[1], _fp);
253 32551 : U1[THMVACE1D::RHOUA] *= dir_mult;
254 :
255 32551 : _fluxes[0] = _numerical_flux.getFlux(_local_side_ids[0], _elem_ids[0], true, U0, U1, _normal[0]);
256 :
257 32551 : _fluxes[1] = _numerical_flux.getFlux(_local_side_ids[0], _elem_ids[0], false, U0, U1, _normal[0]);
258 32551 : _fluxes[1][THMVACE1D::RHOA] *= dir_mult;
259 32551 : _fluxes[1][THMVACE1D::RHOEA] *= dir_mult;
260 97653 : }
261 :
262 : const std::vector<ADReal> &
263 167155 : ADJunctionOneToOne1PhaseUserObject::getFlux(const unsigned int & connection_index) const
264 : {
265 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
266 :
267 167155 : return _fluxes[connection_index];
268 : }
269 :
270 : std::vector<ADReal>
271 56086 : ADJunctionOneToOne1PhaseUserObject::computeElementPrimitiveVariables(const Elem * elem) const
272 : {
273 : const auto U =
274 56086 : FlowModel1PhaseUtils::getElementalSolutionVector<true>(elem, _U_vars, _is_implicit);
275 112172 : return FlowModel1PhaseUtils::computePrimitiveSolutionVector<true>(U, _fp);
276 : }
|
