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 274 : QuadSubChannel1PhaseProblem::validParams()
19 : {
20 274 : InputParameters params = SubChannel1PhaseProblem::validParams();
21 274 : params.addClassDescription(
22 : "Solver class for subchannels in a square lattice assembly and bare fuel pins");
23 274 : return params;
24 0 : }
25 :
26 137 : QuadSubChannel1PhaseProblem::QuadSubChannel1PhaseProblem(const InputParameters & params)
27 137 : : SubChannel1PhaseProblem(params), _subchannel_mesh(SCM::getMesh<QuadSubChannelMesh>(_mesh))
28 : {
29 137 : }
30 :
31 : void
32 196 : QuadSubChannel1PhaseProblem::initializeSolution()
33 : {
34 196 : detectDeformation();
35 :
36 : /// update surface area, wetted perimeter based on: Dpin, displacement
37 196 : 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.setGapWidth(
167 240 : iz, i_gap, (pitch - (*_Dpin_soln)(pin_node_1)) / 2.0 + side_gap + displacement);
168 : }
169 : else // center gap
170 : {
171 120 : _subchannel_mesh.setGapWidth(
172 120 : iz, i_gap, 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.getGapWidth(iz, i_gap) <= 0.0)
176 0 : _subchannel_mesh.setGapWidth(iz, i_gap, 0.0);
177 : }
178 : }
179 5 : }
180 :
181 4813 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
182 : {
183 77354 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
184 : {
185 72737 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
186 72737 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
187 72737 : _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 196 : _aux->solution().close();
194 196 : }
195 :
196 : Real
197 191130 : QuadSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz) const
198 : {
199 191130 : if (!_pin_mesh_exist)
200 : return 0.0;
201 :
202 : // Compute axial location of nodes.
203 185730 : auto z2 = _z_grid[iz];
204 185730 : auto z1 = _z_grid[iz - 1];
205 185730 : auto heated_length = _subchannel_mesh.getHeatedLength();
206 185730 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
207 185730 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
208 183885 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
209 : {
210 : // Compute the height of this element.
211 182040 : auto dz = z2 - z1;
212 : auto heat_rate_in = 0.0;
213 : auto heat_rate_out = 0.0;
214 613960 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
215 : {
216 431920 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
217 431920 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
218 431920 : heat_rate_out += 0.25 * (*_q_prime_soln)(node_out);
219 431920 : heat_rate_in += 0.25 * (*_q_prime_soln)(node_in);
220 : }
221 182040 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
222 : }
223 : else
224 : return 0.0;
225 : }
226 :
227 : Real
228 0 : QuadSubChannel1PhaseProblem::getSubChannelPeripheralDuctWidth(unsigned int i_ch) const
229 : {
230 0 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
231 0 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
232 : {
233 0 : auto width = _subchannel_mesh.getPitch();
234 0 : if (subch_type == EChannelType::CORNER)
235 0 : width += 2.0 * _subchannel_mesh.getSideGap();
236 0 : return width;
237 : }
238 : else
239 0 : mooseError("Channel is not a perimetric subchannel ");
240 : }
241 :
242 : void
243 402 : QuadSubChannel1PhaseProblem::computeh(int iblock)
244 : {
245 402 : unsigned int last_node = (iblock + 1) * _block_size;
246 402 : unsigned int first_node = iblock * _block_size + 1;
247 402 : if (iblock == 0)
248 : {
249 7842 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
250 : {
251 7440 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
252 7440 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
253 7440 : if (h_out < 0)
254 : {
255 0 : mooseError(
256 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
257 : }
258 7440 : _h_soln->set(node, h_out);
259 : }
260 : }
261 :
262 402 : if (!_implicit_bool)
263 : {
264 4994 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
265 : {
266 4810 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
267 34300 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
268 : {
269 29490 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
270 29490 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
271 29490 : auto mdot_in = (*_mdot_soln)(node_in);
272 29490 : auto h_in = (*_h_soln)(node_in); // J/kg
273 29490 : auto volume = dz * (*_S_flow_soln)(node_in);
274 29490 : auto mdot_out = (*_mdot_soln)(node_out);
275 29490 : auto h_out = 0.0;
276 : Real sumWijh = 0.0;
277 : Real sumWijPrimeDhij = 0.0;
278 29490 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
279 : // Calculate sum of crossflow into channel i from channels j around i
280 : unsigned int counter = 0;
281 105690 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
282 : {
283 76200 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
284 : unsigned int ii_ch = chans.first;
285 : // i is always the smallest and first index in the mapping
286 : unsigned int jj_ch = chans.second;
287 76200 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
288 76200 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
289 : // Define donor enthalpy
290 : Real h_star = 0.0;
291 76200 : if (_Wij(i_gap, iz) > 0.0)
292 28008 : h_star = (*_h_soln)(node_in_i);
293 48192 : else if (_Wij(i_gap, iz) < 0.0)
294 24912 : h_star = (*_h_soln)(node_in_j);
295 : // take care of the sign by applying the map, use donor cell
296 76200 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
297 76200 : sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
298 76200 : (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
299 76200 : counter++;
300 : }
301 88470 : h_out = (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy +
302 29490 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
303 29490 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
304 29490 : if (h_out < 0)
305 : {
306 0 : mooseError(name(),
307 : " : Calculation of negative Enthalpy h_out = : ",
308 : h_out,
309 : " Axial Level= : ",
310 : iz);
311 : }
312 29490 : _h_soln->set(node_out, h_out); // J/kg
313 : }
314 : }
315 : }
316 : else
317 : {
318 218 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
319 218 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
320 218 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
321 218 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
322 218 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
323 218 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
324 218 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
325 218 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
326 218 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
327 7368 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
328 : {
329 7150 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
330 7150 : auto iz_ind = iz - first_node;
331 167620 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
332 : {
333 160470 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
334 160470 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
335 160470 : auto S_in = (*_S_flow_soln)(node_in);
336 160470 : auto S_out = (*_S_flow_soln)(node_out);
337 160470 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
338 160470 : auto volume = dz * S_interp;
339 :
340 : // interpolation weight coefficient
341 160470 : auto alpha = computeInterpolationCoefficients(0.5);
342 :
343 : /// Time derivative term
344 160470 : if (iz == first_node)
345 : {
346 : PetscScalar value_vec_tt =
347 5238 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
348 5238 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
349 5238 : LibmeshPetscCall(
350 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
351 : }
352 : else
353 : {
354 155232 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
355 155232 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
356 155232 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
357 155232 : LibmeshPetscCall(MatSetValues(
358 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
359 : }
360 :
361 : // Adding diagonal elements
362 160470 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
363 160470 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
364 160470 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
365 160470 : LibmeshPetscCall(MatSetValues(
366 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
367 :
368 : // Adding RHS elements
369 : PetscScalar rho_old_interp =
370 160470 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), 0.5);
371 : PetscScalar h_old_interp =
372 160470 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), 0.5);
373 160470 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
374 160470 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
375 160470 : LibmeshPetscCall(
376 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
377 :
378 : /// Advective derivative term
379 160470 : if (iz == first_node)
380 : {
381 5238 : PetscScalar value_vec_at = (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
382 5238 : PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
383 5238 : LibmeshPetscCall(VecSetValues(
384 : _hc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
385 : }
386 : else
387 : {
388 155232 : PetscInt row_at = i_ch + _n_channels * iz_ind;
389 155232 : PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
390 155232 : PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in);
391 155232 : LibmeshPetscCall(MatSetValues(
392 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
393 : }
394 :
395 : // Adding diagonal elements
396 160470 : PetscInt row_at = i_ch + _n_channels * iz_ind;
397 160470 : PetscInt col_at = i_ch + _n_channels * iz_ind;
398 160470 : PetscScalar value_at = (*_mdot_soln)(node_out);
399 160470 : LibmeshPetscCall(MatSetValues(
400 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
401 :
402 : /// Cross derivative term
403 : unsigned int counter = 0;
404 : unsigned int cross_index = iz; // iz-1;
405 658470 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
406 : {
407 498000 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
408 : unsigned int ii_ch = chans.first;
409 : unsigned int jj_ch = chans.second;
410 498000 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
411 498000 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
412 : PetscScalar h_star;
413 : // figure out donor axial velocity
414 498000 : if (_Wij(i_gap, cross_index) > 0.0)
415 : {
416 313174 : if (iz == first_node)
417 : {
418 9124 : h_star = (*_h_soln)(node_in_i);
419 18248 : PetscScalar value_vec_ct = -1.0 * alpha *
420 9124 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
421 9124 : _Wij(i_gap, cross_index) * h_star;
422 9124 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
423 9124 : LibmeshPetscCall(VecSetValues(
424 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
425 : }
426 : else
427 : {
428 304050 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
429 304050 : _Wij(i_gap, cross_index);
430 304050 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
431 304050 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
432 304050 : LibmeshPetscCall(MatSetValues(
433 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
434 : }
435 313174 : PetscScalar value_ct = (1.0 - alpha) *
436 313174 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
437 313174 : _Wij(i_gap, cross_index);
438 313174 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
439 313174 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
440 313174 : LibmeshPetscCall(MatSetValues(
441 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
442 : }
443 184826 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
444 : {
445 170426 : if (iz == first_node)
446 : {
447 5996 : h_star = (*_h_soln)(node_in_j);
448 11992 : PetscScalar value_vec_ct = -1.0 * alpha *
449 5996 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
450 5996 : _Wij(i_gap, cross_index) * h_star;
451 5996 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
452 5996 : LibmeshPetscCall(VecSetValues(
453 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
454 : }
455 : else
456 : {
457 164430 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
458 164430 : _Wij(i_gap, cross_index);
459 164430 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
460 164430 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
461 164430 : LibmeshPetscCall(MatSetValues(
462 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
463 : }
464 170426 : PetscScalar value_ct = (1.0 - alpha) *
465 170426 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
466 170426 : _Wij(i_gap, cross_index);
467 170426 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
468 170426 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
469 170426 : LibmeshPetscCall(MatSetValues(
470 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
471 : }
472 :
473 : // Turbulent cross flows
474 498000 : if (iz == first_node)
475 : {
476 : PetscScalar value_vec_ct =
477 16560 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
478 16560 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
479 16560 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
480 16560 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
481 16560 : LibmeshPetscCall(
482 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
483 : }
484 : else
485 : {
486 481440 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
487 481440 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
488 481440 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
489 481440 : LibmeshPetscCall(MatSetValues(
490 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
491 :
492 481440 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
493 481440 : row_ct = i_ch + _n_channels * iz_ind;
494 481440 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
495 481440 : LibmeshPetscCall(MatSetValues(
496 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
497 :
498 481440 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
499 481440 : row_ct = i_ch + _n_channels * iz_ind;
500 481440 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
501 481440 : LibmeshPetscCall(MatSetValues(
502 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
503 : }
504 498000 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
505 498000 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
506 498000 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
507 498000 : LibmeshPetscCall(MatSetValues(
508 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
509 :
510 498000 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
511 498000 : row_ct = i_ch + _n_channels * iz_ind;
512 498000 : col_ct = jj_ch + _n_channels * iz_ind;
513 498000 : LibmeshPetscCall(MatSetValues(
514 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
515 :
516 498000 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
517 498000 : row_ct = i_ch + _n_channels * iz_ind;
518 498000 : col_ct = ii_ch + _n_channels * iz_ind;
519 498000 : LibmeshPetscCall(MatSetValues(
520 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
521 498000 : counter++;
522 : }
523 :
524 : /// Added heat enthalpy
525 160470 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
526 160470 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
527 160470 : LibmeshPetscCall(
528 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
529 : }
530 : }
531 : /// Assembling system
532 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
533 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
534 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
535 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
536 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
537 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
538 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
539 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
540 : // Matrix
541 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
542 218 : LibmeshPetscCall(MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
543 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
544 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
545 218 : LibmeshPetscCall(
546 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
547 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
548 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
549 218 : LibmeshPetscCall(
550 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
551 : #else
552 : LibmeshPetscCall(
553 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
554 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
555 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
556 : LibmeshPetscCall(
557 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
558 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
559 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
560 : LibmeshPetscCall(
561 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
562 : #endif
563 218 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
564 218 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
565 218 : if (_verbose_subchannel)
566 98 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
567 : // RHS
568 218 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
569 218 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
570 218 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
571 218 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
572 :
573 : // Use system to solve for and populate enthalpy
574 218 : LibmeshPetscCall(this->solveAndPopulateEnthalpy(
575 : _hc_sys_h_mat, _hc_sys_h_rhs, first_node, last_node, "h_sys_"));
576 : }
577 402 : }
|
