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 "QuadSubChannel1PhaseProblem.h"
11 : #include "AuxiliarySystem.h"
12 : #include "SCM.h"
13 : #include "SinglePhaseFluidProperties.h"
14 :
15 : registerMooseObject("SubChannelApp", QuadSubChannel1PhaseProblem);
16 :
17 : InputParameters
18 292 : QuadSubChannel1PhaseProblem::validParams()
19 : {
20 292 : InputParameters params = SubChannel1PhaseProblem::validParams();
21 292 : params.addClassDescription(
22 : "Solver class for subchannels in a square lattice assembly and bare fuel pins");
23 292 : return params;
24 0 : }
25 :
26 146 : QuadSubChannel1PhaseProblem::QuadSubChannel1PhaseProblem(const InputParameters & params)
27 146 : : SubChannel1PhaseProblem(params), _subchannel_mesh(SCM::getMesh<QuadSubChannelMesh>(_mesh))
28 : {
29 146 : }
30 :
31 : void
32 200 : QuadSubChannel1PhaseProblem::initializeSolution()
33 : {
34 200 : detectDeformation();
35 :
36 : /// update surface area, wetted perimeter based on: Dpin, displacement
37 200 : if (_deformation)
38 : {
39 5 : _console << " =========== DEFORMATION RECALCULATION ACTIVATED ============== " << std::endl;
40 : Real standard_area, additional_area, wetted_perimeter, displaced_area;
41 5 : auto pitch = _subchannel_mesh.getPitch();
42 5 : auto side_gap = _subchannel_mesh.getSideGap();
43 5 : auto z_blockage = _subchannel_mesh.getZBlockage();
44 5 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
45 5 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
46 35 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
47 : {
48 300 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
49 : {
50 270 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
51 270 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
52 270 : auto Z = _z_grid[iz];
53 : Real rod_area = 0.0;
54 : Real rod_perimeter = 0.0;
55 750 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
56 : {
57 480 : auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
58 480 : rod_area += 0.25 * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
59 480 : rod_perimeter += 0.25 * M_PI * (*_Dpin_soln)(pin_node);
60 : }
61 :
62 270 : if (subch_type == EChannelType::CORNER)
63 : {
64 120 : standard_area = 0.25 * pitch * pitch;
65 120 : displaced_area = (2 * side_gap + pitch) * (*_displacement_soln)(node) / sqrt(2) +
66 120 : (*_displacement_soln)(node) * (*_displacement_soln)(node) / 2;
67 120 : additional_area = pitch * side_gap + side_gap * side_gap;
68 120 : wetted_perimeter =
69 120 : rod_perimeter + pitch + 2 * side_gap + 2 * (*_displacement_soln)(node) / sqrt(2);
70 : }
71 150 : else if (subch_type == EChannelType::EDGE)
72 : {
73 120 : standard_area = 0.5 * pitch * pitch;
74 120 : additional_area = pitch * side_gap;
75 120 : displaced_area = pitch * (*_displacement_soln)(node);
76 120 : wetted_perimeter = rod_perimeter + pitch;
77 : }
78 : else
79 : {
80 30 : standard_area = pitch * pitch;
81 : displaced_area = 0.0;
82 : additional_area = 0.0;
83 : wetted_perimeter = rod_perimeter;
84 : }
85 :
86 : /// Calculate subchannel area
87 270 : auto subchannel_area = displaced_area + standard_area + additional_area - rod_area;
88 :
89 : /// Correct subchannel area and wetted perimeter in case of overlapping pins
90 : auto overlapping_pin_area = 0.0;
91 : auto overlapping_wetted_perimeter = 0.0;
92 990 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
93 : {
94 720 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
95 : auto pin_1 = gap_pins.first;
96 : auto pin_2 = gap_pins.second;
97 720 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
98 720 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
99 720 : auto Diameter1 = (*_Dpin_soln)(pin_node_1);
100 720 : auto Radius1 = Diameter1 / 2.0;
101 720 : auto Diameter2 = (*_Dpin_soln)(pin_node_2);
102 720 : auto Radius2 = Diameter2 / 2.0;
103 720 : auto pitch = _subchannel_mesh.getPitch();
104 :
105 720 : if (pitch < (Radius1 + Radius2)) // overlapping pins
106 : {
107 0 : mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
108 0 : auto cos1 =
109 0 : (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
110 0 : auto cos2 =
111 0 : (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
112 0 : auto angle1 = 2.0 * acos(cos1);
113 0 : auto angle2 = 2.0 * acos(cos2);
114 : // half of the intersecting arc-length
115 0 : overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
116 : // Half of the overlapping area
117 0 : overlapping_pin_area +=
118 0 : 0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
119 0 : 0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
120 0 : (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
121 : }
122 : }
123 270 : subchannel_area += overlapping_pin_area; // correct surface area
124 270 : wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
125 :
126 : /// Apply area reduction on subchannels affected by blockage
127 : auto index = 0;
128 540 : for (const auto & i_blockage : index_blockage)
129 : {
130 270 : if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
131 : {
132 5 : subchannel_area *= reduction_blockage[index];
133 : }
134 270 : index++;
135 : }
136 :
137 270 : _S_flow_soln->set(node, subchannel_area);
138 270 : _w_perim_soln->set(node, wetted_perimeter);
139 : }
140 : }
141 : /// update map of gap between pins (gij) based on: Dpin, displacement
142 35 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
143 : {
144 390 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
145 : {
146 360 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
147 : auto pin_1 = gap_pins.first;
148 : auto pin_2 = gap_pins.second;
149 360 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
150 360 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
151 360 : if (pin_1 == pin_2) // Corner or edge gap
152 : {
153 : auto displacement = 0.0;
154 : auto counter = 0.0;
155 1200 : for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
156 : {
157 960 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
158 960 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
159 960 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
160 : {
161 720 : displacement += (*_displacement_soln)(node);
162 720 : counter += 1.0;
163 : }
164 : }
165 240 : displacement = displacement / counter;
166 240 : _subchannel_mesh._gij_map[iz][i_gap] =
167 240 : (pitch - (*_Dpin_soln)(pin_node_1)) / 2.0 + side_gap + displacement;
168 : }
169 : else // center gap
170 : {
171 120 : _subchannel_mesh._gij_map[iz][i_gap] =
172 120 : pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*_Dpin_soln)(pin_node_2) / 2.0;
173 : }
174 : // if pins come in contact, the gap is zero
175 360 : if (_subchannel_mesh._gij_map[iz][i_gap] <= 0.0)
176 0 : _subchannel_mesh._gij_map[iz][i_gap] = 0.0;
177 : }
178 : }
179 5 : }
180 :
181 4887 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
182 : {
183 78054 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
184 : {
185 73367 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
186 73367 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
187 73367 : _mdot_soln->set(node_out, (*_mdot_soln)(node_in));
188 : }
189 : }
190 :
191 : // We must do a global assembly to make sure data is parallel consistent before we do things
192 : // like compute L2 norms
193 200 : _aux->solution().close();
194 200 : }
195 :
196 : Real
197 195000 : QuadSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz) const
198 : {
199 : // Compute axial location of nodes.
200 195000 : auto z2 = _z_grid[iz];
201 195000 : auto z1 = _z_grid[iz - 1];
202 195000 : auto heated_length = _subchannel_mesh.getHeatedLength();
203 195000 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
204 195000 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
205 191715 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
206 : {
207 : // Compute the height of this element.
208 188430 : auto dz = z2 - z1;
209 188430 : if (_pin_mesh_exist)
210 : {
211 : auto heat_rate_in = 0.0;
212 : auto heat_rate_out = 0.0;
213 601960 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
214 : {
215 424240 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
216 424240 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
217 424240 : heat_rate_out += 0.25 * (*_q_prime_soln)(node_out);
218 424240 : heat_rate_in += 0.25 * (*_q_prime_soln)(node_in);
219 : }
220 177720 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
221 : }
222 : else
223 : {
224 10710 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
225 10710 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
226 10710 : return ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
227 : }
228 : }
229 : else
230 : return 0.0;
231 : }
232 :
233 : Real
234 0 : QuadSubChannel1PhaseProblem::getSubChannelPeripheralDuctWidth(unsigned int i_ch) const
235 : {
236 0 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
237 0 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
238 : {
239 0 : auto width = _subchannel_mesh.getPitch();
240 0 : if (subch_type == EChannelType::CORNER)
241 0 : width += 2.0 * _subchannel_mesh.getSideGap();
242 0 : return width;
243 : }
244 : else
245 0 : mooseError("Channel is not a perimetric subchannel ");
246 : }
247 :
248 : void
249 423 : QuadSubChannel1PhaseProblem::computeh(int iblock)
250 : {
251 423 : unsigned int last_node = (iblock + 1) * _block_size;
252 423 : unsigned int first_node = iblock * _block_size + 1;
253 423 : if (iblock == 0)
254 : {
255 8052 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
256 : {
257 7629 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
258 7629 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
259 7629 : if (h_out < 0)
260 : {
261 0 : mooseError(
262 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
263 : }
264 7629 : _h_soln->set(node, h_out);
265 : }
266 : }
267 :
268 423 : if (!_implicit_bool)
269 : {
270 4895 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
271 : {
272 4720 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
273 33400 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
274 : {
275 28680 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
276 28680 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
277 28680 : auto mdot_in = (*_mdot_soln)(node_in);
278 28680 : auto h_in = (*_h_soln)(node_in); // J/kg
279 28680 : auto volume = dz * (*_S_flow_soln)(node_in);
280 28680 : auto mdot_out = (*_mdot_soln)(node_out);
281 28680 : auto h_out = 0.0;
282 : Real sumWijh = 0.0;
283 : Real sumWijPrimeDhij = 0.0;
284 28680 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
285 : // Calculate sum of crossflow into channel i from channels j around i
286 : unsigned int counter = 0;
287 102720 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
288 : {
289 74040 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
290 : unsigned int ii_ch = chans.first;
291 : // i is always the smallest and first index in the mapping
292 : unsigned int jj_ch = chans.second;
293 74040 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
294 74040 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
295 : // Define donor enthalpy
296 : Real h_star = 0.0;
297 74040 : if (_Wij(i_gap, iz) > 0.0)
298 28008 : h_star = (*_h_soln)(node_in_i);
299 46032 : else if (_Wij(i_gap, iz) < 0.0)
300 24912 : h_star = (*_h_soln)(node_in_j);
301 : // take care of the sign by applying the map, use donor cell
302 74040 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
303 74040 : sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
304 74040 : (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
305 74040 : counter++;
306 : }
307 86040 : h_out = (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy +
308 28680 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
309 28680 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
310 28680 : if (h_out < 0)
311 : {
312 0 : mooseError(name(),
313 : " : Calculation of negative Enthalpy h_out = : ",
314 : h_out,
315 : " Axial Level= : ",
316 : iz);
317 : }
318 28680 : _h_soln->set(node_out, h_out); // J/kg
319 : }
320 : }
321 : }
322 : else
323 : {
324 248 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
325 248 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
326 248 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
327 248 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
328 248 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
329 248 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
330 248 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
331 248 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
332 248 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
333 7788 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
334 : {
335 7540 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
336 7540 : auto iz_ind = iz - first_node;
337 171520 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
338 : {
339 163980 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
340 163980 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
341 163980 : auto S_in = (*_S_flow_soln)(node_in);
342 163980 : auto S_out = (*_S_flow_soln)(node_out);
343 163980 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
344 163980 : auto volume = dz * S_interp;
345 :
346 : // interpolation weight coefficient
347 163980 : auto alpha = computeInterpolationCoefficients(0.5);
348 :
349 : /// Time derivative term
350 163980 : if (iz == first_node)
351 : {
352 : PetscScalar value_vec_tt =
353 5508 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
354 5508 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
355 5508 : LibmeshPetscCall(
356 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
357 : }
358 : else
359 : {
360 158472 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
361 158472 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
362 158472 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
363 158472 : LibmeshPetscCall(MatSetValues(
364 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
365 : }
366 :
367 : // Adding diagonal elements
368 163980 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
369 163980 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
370 163980 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
371 163980 : LibmeshPetscCall(MatSetValues(
372 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
373 :
374 : // Adding RHS elements
375 : PetscScalar rho_old_interp =
376 163980 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), 0.5);
377 : PetscScalar h_old_interp =
378 163980 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), 0.5);
379 163980 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
380 163980 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
381 163980 : LibmeshPetscCall(
382 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
383 :
384 : /// Advective derivative term
385 163980 : if (iz == first_node)
386 : {
387 5508 : PetscScalar value_vec_at = (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
388 5508 : PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
389 5508 : LibmeshPetscCall(VecSetValues(
390 : _hc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
391 : }
392 : else
393 : {
394 158472 : PetscInt row_at = i_ch + _n_channels * iz_ind;
395 158472 : PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
396 158472 : PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in);
397 158472 : LibmeshPetscCall(MatSetValues(
398 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
399 : }
400 :
401 : // Adding diagonal elements
402 163980 : PetscInt row_at = i_ch + _n_channels * iz_ind;
403 163980 : PetscInt col_at = i_ch + _n_channels * iz_ind;
404 163980 : PetscScalar value_at = (*_mdot_soln)(node_out);
405 163980 : LibmeshPetscCall(MatSetValues(
406 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
407 :
408 : /// Cross derivative term
409 : unsigned int counter = 0;
410 : unsigned int cross_index = iz; // iz-1;
411 671340 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
412 : {
413 507360 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
414 : unsigned int ii_ch = chans.first;
415 : unsigned int jj_ch = chans.second;
416 507360 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
417 507360 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
418 : PetscScalar h_star;
419 : // figure out donor axial velocity
420 507360 : if (_Wij(i_gap, cross_index) > 0.0)
421 : {
422 317864 : if (iz == first_node)
423 : {
424 9524 : h_star = (*_h_soln)(node_in_i);
425 19048 : PetscScalar value_vec_ct = -1.0 * alpha *
426 9524 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
427 9524 : _Wij(i_gap, cross_index) * h_star;
428 9524 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
429 9524 : LibmeshPetscCall(VecSetValues(
430 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
431 : }
432 : else
433 : {
434 308340 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
435 308340 : _Wij(i_gap, cross_index);
436 308340 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
437 308340 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
438 308340 : LibmeshPetscCall(MatSetValues(
439 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
440 : }
441 317864 : PetscScalar value_ct = (1.0 - alpha) *
442 317864 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
443 317864 : _Wij(i_gap, cross_index);
444 317864 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
445 317864 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
446 317864 : LibmeshPetscCall(MatSetValues(
447 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
448 : }
449 189496 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
450 : {
451 175096 : if (iz == first_node)
452 : {
453 6316 : h_star = (*_h_soln)(node_in_j);
454 12632 : PetscScalar value_vec_ct = -1.0 * alpha *
455 6316 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
456 6316 : _Wij(i_gap, cross_index) * h_star;
457 6316 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
458 6316 : LibmeshPetscCall(VecSetValues(
459 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
460 : }
461 : else
462 : {
463 168780 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
464 168780 : _Wij(i_gap, cross_index);
465 168780 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
466 168780 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
467 168780 : LibmeshPetscCall(MatSetValues(
468 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
469 : }
470 175096 : PetscScalar value_ct = (1.0 - alpha) *
471 175096 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
472 175096 : _Wij(i_gap, cross_index);
473 175096 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
474 175096 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
475 175096 : LibmeshPetscCall(MatSetValues(
476 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
477 : }
478 :
479 : // Turbulent cross flows
480 507360 : if (iz == first_node)
481 : {
482 : PetscScalar value_vec_ct =
483 17280 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
484 17280 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
485 17280 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
486 17280 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
487 17280 : LibmeshPetscCall(
488 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
489 : }
490 : else
491 : {
492 490080 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
493 490080 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
494 490080 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
495 490080 : LibmeshPetscCall(MatSetValues(
496 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
497 :
498 490080 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
499 490080 : row_ct = i_ch + _n_channels * iz_ind;
500 490080 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
501 490080 : LibmeshPetscCall(MatSetValues(
502 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
503 :
504 490080 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
505 490080 : row_ct = i_ch + _n_channels * iz_ind;
506 490080 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
507 490080 : LibmeshPetscCall(MatSetValues(
508 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
509 : }
510 507360 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
511 507360 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
512 507360 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
513 507360 : LibmeshPetscCall(MatSetValues(
514 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
515 :
516 507360 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
517 507360 : row_ct = i_ch + _n_channels * iz_ind;
518 507360 : col_ct = jj_ch + _n_channels * iz_ind;
519 507360 : LibmeshPetscCall(MatSetValues(
520 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
521 :
522 507360 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
523 507360 : row_ct = i_ch + _n_channels * iz_ind;
524 507360 : col_ct = ii_ch + _n_channels * iz_ind;
525 507360 : LibmeshPetscCall(MatSetValues(
526 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
527 507360 : counter++;
528 : }
529 :
530 : /// Added heat enthalpy
531 163980 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
532 163980 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
533 163980 : LibmeshPetscCall(
534 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
535 : }
536 : }
537 : /// Assembling system
538 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
539 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
540 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
541 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
542 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
543 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
544 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
545 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
546 : // Matrix
547 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
548 248 : LibmeshPetscCall(MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
549 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
550 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
551 248 : LibmeshPetscCall(
552 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
553 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
554 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
555 248 : LibmeshPetscCall(
556 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
557 : #else
558 : LibmeshPetscCall(
559 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
560 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
561 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
562 : LibmeshPetscCall(
563 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
564 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
565 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
566 : LibmeshPetscCall(
567 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
568 : #endif
569 248 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
570 248 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
571 248 : if (_verbose_subchannel)
572 128 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
573 : // RHS
574 248 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
575 248 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
576 248 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
577 248 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
578 :
579 : // Use system to solve for and populate enthalpy
580 248 : LibmeshPetscCall(this->solveAndPopulateEnthalpy(
581 : _hc_sys_h_mat, _hc_sys_h_rhs, first_node, last_node, "h_sys_"));
582 : }
583 423 : }
|
