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 :
14 : registerMooseObject("SubChannelApp", QuadSubChannel1PhaseProblem);
15 :
16 : InputParameters
17 284 : QuadSubChannel1PhaseProblem::validParams()
18 : {
19 284 : InputParameters params = SubChannel1PhaseProblem::validParams();
20 284 : params.addClassDescription(
21 : "Solver class for subchannels in a square lattice assembly and bare fuel pins");
22 568 : params.addRequiredParam<Real>("beta",
23 : "Thermal diffusion coefficient used in turbulent crossflow.");
24 568 : params.addParam<bool>(
25 : "default_friction_model",
26 568 : true,
27 : "Boolean to define which friction model to use (default: Pang, B. et al. "
28 : "KIT, 2013. / non-default: Todreas-Kazimi NUCLEAR SYSTEMS, second edition, Volume 1, 2011)");
29 568 : params.addParam<bool>(
30 : "constant_beta",
31 568 : true,
32 : "Boolean to define the use of a constant beta or beta correlation (Kim and Chung, 2001)");
33 284 : return params;
34 0 : }
35 :
36 142 : QuadSubChannel1PhaseProblem::QuadSubChannel1PhaseProblem(const InputParameters & params)
37 : : SubChannel1PhaseProblem(params),
38 284 : _subchannel_mesh(SCM::getMesh<QuadSubChannelMesh>(_mesh)),
39 284 : _beta(getParam<Real>("beta")),
40 284 : _default_friction_model(getParam<bool>("default_friction_model")),
41 426 : _constant_beta(getParam<bool>("constant_beta"))
42 : {
43 142 : }
44 :
45 : void
46 190 : QuadSubChannel1PhaseProblem::initializeSolution()
47 : {
48 : /// update surface area, wetted perimeter based on: Dpin, displacement
49 190 : if (_deformation)
50 : {
51 : Real standard_area, additional_area, wetted_perimeter, displaced_area;
52 6 : auto pitch = _subchannel_mesh.getPitch();
53 6 : auto side_gap = _subchannel_mesh.getSideGap();
54 6 : auto z_blockage = _subchannel_mesh.getZBlockage();
55 6 : auto index_blockage = _subchannel_mesh.getIndexBlockage();
56 6 : auto reduction_blockage = _subchannel_mesh.getReductionBlockage();
57 132 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
58 : {
59 4662 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
60 : {
61 4536 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
62 4536 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
63 4536 : auto Z = _z_grid[iz];
64 : Real rod_area = 0.0;
65 : Real rod_perimeter = 0.0;
66 17136 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
67 : {
68 12600 : auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
69 12600 : rod_area += 0.25 * 0.25 * M_PI * (*_Dpin_soln)(pin_node) * (*_Dpin_soln)(pin_node);
70 12600 : rod_perimeter += 0.25 * M_PI * (*_Dpin_soln)(pin_node);
71 : }
72 :
73 4536 : if (subch_type == EChannelType::CORNER)
74 : {
75 504 : standard_area = 0.25 * pitch * pitch;
76 504 : displaced_area = (2 * side_gap + pitch) * (*_displacement_soln)(node) / sqrt(2) +
77 504 : (*_displacement_soln)(node) * (*_displacement_soln)(node) / 2;
78 504 : additional_area = pitch * side_gap + side_gap * side_gap;
79 504 : wetted_perimeter =
80 504 : rod_perimeter + pitch + 2 * side_gap + 2 * (*_displacement_soln)(node) / sqrt(2);
81 : }
82 4032 : else if (subch_type == EChannelType::EDGE)
83 : {
84 2016 : standard_area = 0.5 * pitch * pitch;
85 2016 : additional_area = pitch * side_gap;
86 2016 : displaced_area = pitch * (*_displacement_soln)(node);
87 2016 : wetted_perimeter = rod_perimeter + pitch;
88 : }
89 : else
90 : {
91 2016 : standard_area = pitch * pitch;
92 : displaced_area = 0.0;
93 : additional_area = 0.0;
94 : wetted_perimeter = rod_perimeter;
95 : }
96 :
97 : /// Calculate subchannel area
98 4536 : auto subchannel_area = displaced_area + standard_area + additional_area - rod_area;
99 :
100 : /// Correct subchannel area and wetted perimeter in case of overlapping pins
101 : auto overlapping_pin_area = 0.0;
102 : auto overlapping_wetted_perimeter = 0.0;
103 19656 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
104 : {
105 15120 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
106 : auto pin_1 = gap_pins.first;
107 : auto pin_2 = gap_pins.second;
108 15120 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
109 15120 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
110 15120 : auto Diameter1 = (*_Dpin_soln)(pin_node_1);
111 15120 : auto Radius1 = Diameter1 / 2.0;
112 15120 : auto Diameter2 = (*_Dpin_soln)(pin_node_2);
113 15120 : auto Radius2 = Diameter2 / 2.0;
114 15120 : auto pitch = _subchannel_mesh.getPitch();
115 :
116 15120 : if (pitch < (Radius1 + Radius2)) // overlapping pins
117 : {
118 0 : mooseWarning(" The gap of index : '", i_gap, " at axial cell ", iz, " ' is blocked.");
119 0 : auto cos1 =
120 0 : (pitch * pitch + Radius1 * Radius1 - Radius2 * Radius2) / (2 * pitch * Radius1);
121 0 : auto cos2 =
122 0 : (pitch * pitch + Radius2 * Radius2 - Radius1 * Radius1) / (2 * pitch * Radius2);
123 0 : auto angle1 = 2.0 * acos(cos1);
124 0 : auto angle2 = 2.0 * acos(cos2);
125 : // half of the intersecting arc-length
126 0 : overlapping_wetted_perimeter += 0.5 * angle1 * Radius1 + 0.5 * angle2 * Radius2;
127 : // Half of the overlapping area
128 0 : overlapping_pin_area +=
129 0 : 0.5 * Radius1 * Radius1 * acos(cos1) + 0.5 * Radius2 * Radius2 * acos(cos2) -
130 0 : 0.25 * sqrt((-pitch + Radius1 + Radius2) * (pitch + Radius1 - Radius2) *
131 0 : (pitch - Radius1 + Radius2) * (pitch + Radius1 + Radius2));
132 : }
133 : }
134 4536 : subchannel_area += overlapping_pin_area; // correct surface area
135 4536 : wetted_perimeter += -overlapping_wetted_perimeter; // correct wetted perimeter
136 :
137 : /// Apply area reduction on subchannels affected by blockage
138 : auto index = 0;
139 9072 : for (const auto & i_blockage : index_blockage)
140 : {
141 4536 : if (i_ch == i_blockage && (Z >= z_blockage.front() && Z <= z_blockage.back()))
142 : {
143 6 : subchannel_area *= reduction_blockage[index];
144 : }
145 4536 : index++;
146 : }
147 :
148 4536 : _S_flow_soln->set(node, subchannel_area);
149 4536 : _w_perim_soln->set(node, wetted_perimeter);
150 : }
151 : }
152 : /// update map of gap between pins (gij) based on: Dpin, displacement
153 132 : for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
154 : {
155 7686 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
156 : {
157 7560 : auto gap_pins = _subchannel_mesh.getGapPins(i_gap);
158 : auto pin_1 = gap_pins.first;
159 : auto pin_2 = gap_pins.second;
160 7560 : auto * pin_node_1 = _subchannel_mesh.getPinNode(pin_1, iz);
161 7560 : auto * pin_node_2 = _subchannel_mesh.getPinNode(pin_2, iz);
162 7560 : if (pin_1 == pin_2) // Corner or edge gap
163 : {
164 : auto displacement = 0.0;
165 : auto counter = 0.0;
166 12600 : for (auto i_ch : _subchannel_mesh.getPinChannels(pin_1))
167 : {
168 10080 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
169 10080 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
170 10080 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
171 : {
172 6048 : displacement += (*_displacement_soln)(node);
173 6048 : counter += 1.0;
174 : }
175 : }
176 2520 : displacement = displacement / counter;
177 2520 : _subchannel_mesh._gij_map[iz][i_gap] =
178 2520 : (pitch - (*_Dpin_soln)(pin_node_1)) / 2.0 + side_gap + displacement;
179 : }
180 : else // center gap
181 : {
182 5040 : _subchannel_mesh._gij_map[iz][i_gap] =
183 5040 : pitch - (*_Dpin_soln)(pin_node_1) / 2.0 - (*_Dpin_soln)(pin_node_2) / 2.0;
184 : }
185 : // if pins come in contact, the gap is zero
186 7560 : if (_subchannel_mesh._gij_map[iz][i_gap] <= 0.0)
187 0 : _subchannel_mesh._gij_map[iz][i_gap] = 0.0;
188 : }
189 : }
190 6 : }
191 :
192 7180 : for (unsigned int iz = 1; iz < _n_cells + 1; iz++)
193 : {
194 118856 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
195 : {
196 111866 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
197 111866 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
198 111866 : _mdot_soln->set(node_out, (*_mdot_soln)(node_in));
199 : }
200 : }
201 :
202 : // We must do a global assembly to make sure data is parallel consistent before we do things
203 : // like compute L2 norms
204 190 : _aux->solution().close();
205 190 : }
206 :
207 : Real
208 13201524 : QuadSubChannel1PhaseProblem::computeFrictionFactor(FrictionStruct friction_args)
209 : {
210 13201524 : auto Re = friction_args.Re;
211 13201524 : auto i_ch = friction_args.i_ch;
212 : /// Pang, B. et al. KIT, 2013
213 13201524 : if (_default_friction_model)
214 : {
215 : Real a, b;
216 10780404 : if (Re < 1)
217 : {
218 : return 64.0;
219 : }
220 10780404 : else if (Re >= 1 and Re < 5000)
221 : {
222 : a = 64.0;
223 : b = -1.0;
224 : }
225 8072964 : else if (Re >= 5000 and Re < 30000)
226 : {
227 : a = 0.316;
228 : b = -0.25;
229 : }
230 : else
231 : {
232 : a = 0.184;
233 : b = -0.20;
234 : }
235 10780404 : return a * std::pow(Re, b);
236 : }
237 : /// Todreas-Kazimi NUCLEAR SYSTEMS, second edition, Volume 1, 2011
238 : else
239 : {
240 : Real aL, b1L, b2L, cL;
241 : Real aT, b1T, b2T, cT;
242 2421120 : auto pitch = _subchannel_mesh.getPitch();
243 2421120 : auto pin_diameter = _subchannel_mesh.getPinDiameter();
244 : // This gap is a constant value for the whole assembly. Might want to make it
245 : // subchannel specific in the future if we have duct deformation.
246 2421120 : auto side_gap = _subchannel_mesh.getSideGap();
247 2421120 : auto w = (pin_diameter / 2.0) + (pitch / 2.0) + side_gap;
248 2421120 : auto p_over_d = pitch / pin_diameter;
249 2421120 : auto w_over_d = w / pin_diameter;
250 2421120 : auto ReL = std::pow(10, (p_over_d - 1)) * 320.0;
251 2421120 : auto ReT = std::pow(10, 0.7 * (p_over_d - 1)) * 1.0E+4;
252 2421120 : auto psi = std::log(Re / ReL) / std::log(ReT / ReL);
253 2421120 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
254 : const Real lambda = 7.0;
255 :
256 : // Find the coefficients of bare Pin bundle friction factor
257 : // correlations for turbulent and laminar flow regimes. Todreas & Kazimi, Nuclear Systems Volume
258 : // 1
259 2421120 : if (subch_type == EChannelType::CENTER)
260 : {
261 1116720 : if (p_over_d < 1.1)
262 : {
263 : aL = 26.37;
264 : b1L = 374.2;
265 : b2L = -493.9;
266 : aT = 0.09423;
267 : b1T = 0.5806;
268 : b2T = -1.239;
269 : }
270 : else
271 : {
272 : aL = 35.55;
273 : b1L = 263.7;
274 : b2L = -190.2;
275 : aT = 0.1339;
276 : b1T = 0.09059;
277 : b2T = -0.09926;
278 : }
279 : // laminar flow friction factor for bare Pin bundle - Center subchannel
280 1116720 : cL = aL + b1L * (p_over_d - 1) + b2L * std::pow((p_over_d - 1), 2);
281 : // turbulent flow friction factor for bare Pin bundle - Center subchannel
282 1116720 : cT = aT + b1T * (p_over_d - 1) + b2T * std::pow((p_over_d - 1), 2);
283 : }
284 1304400 : else if (subch_type == EChannelType::EDGE)
285 : {
286 1043520 : if (p_over_d < 1.1)
287 : {
288 : aL = 26.18;
289 : b1L = 554.5;
290 : b2L = -1480;
291 : aT = 0.09377;
292 : b1T = 0.8732;
293 : b2T = -3.341;
294 : }
295 : else
296 : {
297 : aL = 44.40;
298 : b1L = 256.7;
299 : b2L = -267.6;
300 : aT = 0.1430;
301 : b1T = 0.04199;
302 : b2T = -0.04428;
303 : }
304 : // laminar flow friction factor for bare Pin bundle - Edge subchannel
305 1043520 : cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2);
306 : // turbulent flow friction factor for bare Pin bundle - Edge subchannel
307 1043520 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
308 : }
309 : else
310 : {
311 260880 : if (p_over_d < 1.1)
312 : {
313 : aL = 28.62;
314 : b1L = 715.9;
315 : b2L = -2807;
316 : aT = 0.09755;
317 : b1T = 1.127;
318 : b2T = -6.304;
319 : }
320 : else
321 : {
322 : aL = 58.83;
323 : b1L = 160.7;
324 : b2L = -203.5;
325 : aT = 0.1452;
326 : b1T = 0.02681;
327 : b2T = -0.03411;
328 : }
329 : // laminar flow friction factor for bare Pin bundle - Corner subchannel
330 260880 : cL = aL + b1L * (w_over_d - 1) + b2L * std::pow((w_over_d - 1), 2);
331 : // turbulent flow friction factor for bare Pin bundle - Corner subchannel
332 260880 : cT = aT + b1T * (w_over_d - 1) + b2T * std::pow((w_over_d - 1), 2);
333 : }
334 :
335 : // laminar friction factor
336 2421120 : auto fL = cL / Re;
337 : // turbulent friction factor
338 2421120 : auto fT = cT / std::pow(Re, 0.18);
339 :
340 2421120 : if (Re < ReL)
341 : {
342 : // laminar flow
343 : return fL;
344 : }
345 2421120 : else if (Re > ReT)
346 : {
347 : // turbulent flow
348 : return fT;
349 : }
350 : else
351 : {
352 : // transient flow: psi definition uses a Bulk ReT/ReL number, same for all channels
353 0 : return fL * std::pow((1 - psi), 1.0 / 3.0) * (1 - std::pow(psi, lambda)) +
354 0 : fT * std::pow(psi, 1.0 / 3.0);
355 : }
356 : }
357 : }
358 :
359 : Real
360 21571344 : QuadSubChannel1PhaseProblem::computeBeta(unsigned int i_gap, unsigned int iz, bool /*enthalpy*/)
361 : {
362 21571344 : auto beta = _beta;
363 21571344 : if (!_constant_beta)
364 : {
365 3913200 : const Real & pitch = _subchannel_mesh.getPitch();
366 3913200 : const Real & pin_diameter = _subchannel_mesh.getPinDiameter();
367 3913200 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
368 : unsigned int i_ch = chans.first;
369 : unsigned int j_ch = chans.second;
370 3913200 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
371 3913200 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
372 3913200 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
373 3913200 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
374 3913200 : auto Si_in = (*_S_flow_soln)(node_in_i);
375 3913200 : auto Sj_in = (*_S_flow_soln)(node_in_j);
376 3913200 : auto Si_out = (*_S_flow_soln)(node_out_i);
377 3913200 : auto Sj_out = (*_S_flow_soln)(node_out_j);
378 : // crossflow area between channels i,j (dz*gap_width)
379 3913200 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
380 : auto avg_massflux =
381 3913200 : 0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
382 3913200 : ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
383 3913200 : auto S_total = Si_in + Sj_in + Si_out + Sj_out;
384 3913200 : auto Si = 0.5 * (Si_in + Si_out);
385 3913200 : auto Sj = 0.5 * (Sj_in + Sj_out);
386 3913200 : auto w_perim_i = 0.5 * ((*_w_perim_soln)(node_in_i) + (*_w_perim_soln)(node_out_i));
387 3913200 : auto w_perim_j = 0.5 * ((*_w_perim_soln)(node_in_j) + (*_w_perim_soln)(node_out_j));
388 3913200 : auto avg_mu = (1 / S_total) * ((*_mu_soln)(node_out_i)*Si_out + (*_mu_soln)(node_in_i)*Si_in +
389 3913200 : (*_mu_soln)(node_out_j)*Sj_out + (*_mu_soln)(node_in_j)*Sj_in);
390 3913200 : auto avg_hD = 4.0 * (Si + Sj) / (w_perim_i + w_perim_j);
391 3913200 : auto Re = avg_massflux * avg_hD / avg_mu;
392 : Real gamma = 20.0; // empirical constant
393 : Real sf = 1.0; // shape factor
394 : Real a = 0.18;
395 : Real b = 0.2;
396 3913200 : auto f = a * std::pow(Re, -b); // Rehme 1992 circular tube friction factor
397 : auto k = (1 / S_total) *
398 3913200 : (_fp->k_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
399 3913200 : _fp->k_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
400 3913200 : _fp->k_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
401 3913200 : _fp->k_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
402 : auto cp = (1 / S_total) *
403 3913200 : (_fp->cp_from_p_T((*_P_soln)(node_out_i) + _P_out, (*_T_soln)(node_out_i)) * Si_out +
404 3913200 : _fp->cp_from_p_T((*_P_soln)(node_in_i) + _P_out, (*_T_soln)(node_in_i)) * Si_in +
405 3913200 : _fp->cp_from_p_T((*_P_soln)(node_out_j) + _P_out, (*_T_soln)(node_out_j)) * Sj_out +
406 3913200 : _fp->cp_from_p_T((*_P_soln)(node_in_j) + _P_out, (*_T_soln)(node_in_j)) * Sj_in);
407 3913200 : auto Pr = avg_mu * cp / k; // Prandtl number
408 3913200 : auto Pr_t = Pr * (Re / gamma) * std::sqrt(f / 8.0); // Turbulent Prandtl number
409 3913200 : auto delta = pitch; // centroid to centroid distance
410 : auto L_x = sf * delta; // axial length scale (gap is the lateral length scale)
411 3913200 : auto lamda = gap / L_x; // aspect ratio
412 3913200 : auto a_x = 1.0 - 2.0 * lamda * lamda / libMesh::pi; // velocity coefficient
413 3913200 : auto z_FP_over_D = (2.0 * L_x / pin_diameter) *
414 3913200 : (1 + (-0.5 * std::log(lamda) + 0.5 * std::log(4.0) - 0.25) * lamda * lamda);
415 3913200 : auto Str = 1.0 / (0.822 * (gap / pin_diameter) + 0.144); // Strouhal number (Wu & Trupp 1994)
416 3913200 : auto freq_factor = 2.0 / std::pow(gamma, 2) * std::sqrt(a / 8.0) * (avg_hD / gap);
417 3913200 : auto rod_mixing = (1 / Pr_t) * lamda;
418 3913200 : auto axial_mixing = a_x * z_FP_over_D * Str;
419 : // Mixing Stanton number: Stg (eq 25,Kim and Chung (2001), eq 19 (Jeong et. al 2005)
420 3913200 : beta = freq_factor * (rod_mixing + axial_mixing) * std::pow(Re, -b / 2.0);
421 : }
422 : mooseAssert(beta > 0, "beta should be positive");
423 21571344 : return beta;
424 : }
425 :
426 : Real
427 328272 : QuadSubChannel1PhaseProblem::computeAddedHeatPin(unsigned int i_ch, unsigned int iz)
428 : {
429 : // Compute axial location of nodes.
430 328272 : auto z2 = _z_grid[iz];
431 328272 : auto z1 = _z_grid[iz - 1];
432 328272 : auto heated_length = _subchannel_mesh.getHeatedLength();
433 328272 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
434 328272 : if (MooseUtils::absoluteFuzzyGreaterThan(z2, unheated_length_entry) &&
435 324006 : MooseUtils::absoluteFuzzyLessThan(z1, unheated_length_entry + heated_length))
436 : {
437 : // Compute the height of this element.
438 319740 : auto dz = z2 - z1;
439 319740 : if (_pin_mesh_exist)
440 : {
441 : auto heat_rate_in = 0.0;
442 : auto heat_rate_out = 0.0;
443 1044720 : for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
444 : {
445 743520 : auto * node_in = _subchannel_mesh.getPinNode(i_pin, iz - 1);
446 743520 : auto * node_out = _subchannel_mesh.getPinNode(i_pin, iz);
447 743520 : heat_rate_out += 0.25 * (*_q_prime_soln)(node_out);
448 743520 : heat_rate_in += 0.25 * (*_q_prime_soln)(node_in);
449 : }
450 301200 : return (heat_rate_in + heat_rate_out) * dz / 2.0;
451 : }
452 : else
453 : {
454 18540 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
455 18540 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
456 18540 : return ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
457 : }
458 : }
459 : else
460 : return 0.0;
461 : }
462 :
463 : void
464 668 : QuadSubChannel1PhaseProblem::computeh(int iblock)
465 : {
466 668 : unsigned int last_node = (iblock + 1) * _block_size;
467 668 : unsigned int first_node = iblock * _block_size + 1;
468 668 : if (iblock == 0)
469 : {
470 15794 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
471 : {
472 15126 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
473 15126 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
474 15126 : if (h_out < 0)
475 : {
476 0 : mooseError(
477 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
478 : }
479 15126 : _h_soln->set(node, h_out);
480 : }
481 : }
482 :
483 668 : if (!_implicit_bool)
484 : {
485 7174 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
486 : {
487 6980 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
488 45380 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
489 : {
490 38400 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
491 38400 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
492 38400 : auto mdot_in = (*_mdot_soln)(node_in);
493 38400 : auto h_in = (*_h_soln)(node_in); // J/kg
494 38400 : auto volume = dz * (*_S_flow_soln)(node_in);
495 38400 : auto mdot_out = (*_mdot_soln)(node_out);
496 38400 : auto h_out = 0.0;
497 : Real sumWijh = 0.0;
498 : Real sumWijPrimeDhij = 0.0;
499 38400 : Real added_enthalpy = computeAddedHeatPin(i_ch, iz);
500 : // Calculate sum of crossflow into channel i from channels j around i
501 : unsigned int counter = 0;
502 135120 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
503 : {
504 96720 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
505 : unsigned int ii_ch = chans.first;
506 : // i is always the smallest and first index in the mapping
507 : unsigned int jj_ch = chans.second;
508 96720 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
509 96720 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
510 : // Define donor enthalpy
511 : Real h_star = 0.0;
512 96720 : if (_Wij(i_gap, iz) > 0.0)
513 42408 : h_star = (*_h_soln)(node_in_i);
514 54312 : else if (_Wij(i_gap, iz) < 0.0)
515 35352 : h_star = (*_h_soln)(node_in_j);
516 : // take care of the sign by applying the map, use donor cell
517 96720 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
518 96720 : sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) -
519 96720 : (*_h_soln)(node_in_j) - (*_h_soln)(node_in_i));
520 96720 : counter++;
521 : }
522 115200 : h_out = (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy +
523 38400 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
524 38400 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
525 38400 : if (h_out < 0)
526 : {
527 0 : mooseError(name(),
528 : " : Calculation of negative Enthalpy h_out = : ",
529 : h_out,
530 : " Axial Level= : ",
531 : iz);
532 : }
533 38400 : _h_soln->set(node_out, h_out); // J/kg
534 : }
535 : }
536 : }
537 : else
538 : {
539 474 : LibmeshPetscCall(MatZeroEntries(_hc_time_derivative_mat));
540 474 : LibmeshPetscCall(MatZeroEntries(_hc_advective_derivative_mat));
541 474 : LibmeshPetscCall(MatZeroEntries(_hc_cross_derivative_mat));
542 474 : LibmeshPetscCall(VecZeroEntries(_hc_time_derivative_rhs));
543 474 : LibmeshPetscCall(VecZeroEntries(_hc_advective_derivative_rhs));
544 474 : LibmeshPetscCall(VecZeroEntries(_hc_cross_derivative_rhs));
545 474 : LibmeshPetscCall(VecZeroEntries(_hc_added_heat_rhs));
546 474 : LibmeshPetscCall(MatZeroEntries(_hc_sys_h_mat));
547 474 : LibmeshPetscCall(VecZeroEntries(_hc_sys_h_rhs));
548 12750 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
549 : {
550 12276 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
551 12276 : auto iz_ind = iz - first_node;
552 299340 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
553 : {
554 287064 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
555 287064 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
556 287064 : auto S_in = (*_S_flow_soln)(node_in);
557 287064 : auto S_out = (*_S_flow_soln)(node_out);
558 287064 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
559 287064 : auto volume = dz * S_interp;
560 :
561 : // interpolation weight coefficient
562 287064 : auto alpha = computeInterpolationCoefficients(0.5);
563 :
564 : /// Time derivative term
565 287064 : if (iz == first_node)
566 : {
567 : PetscScalar value_vec_tt =
568 12420 : -1.0 * _TR * alpha * (*_rho_soln)(node_in) * (*_h_soln)(node_in)*volume / _dt;
569 12420 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
570 12420 : LibmeshPetscCall(
571 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
572 : }
573 : else
574 : {
575 274644 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
576 274644 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
577 274644 : PetscScalar value_tt = _TR * alpha * (*_rho_soln)(node_in)*volume / _dt;
578 274644 : LibmeshPetscCall(MatSetValues(
579 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
580 : }
581 :
582 : // Adding diagonal elements
583 287064 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
584 287064 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
585 287064 : PetscScalar value_tt = _TR * (1.0 - alpha) * (*_rho_soln)(node_out)*volume / _dt;
586 287064 : LibmeshPetscCall(MatSetValues(
587 : _hc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
588 :
589 : // Adding RHS elements
590 : PetscScalar rho_old_interp =
591 287064 : computeInterpolatedValue(_rho_soln->old(node_out), _rho_soln->old(node_in), 0.5);
592 : PetscScalar h_old_interp =
593 287064 : computeInterpolatedValue(_h_soln->old(node_out), _h_soln->old(node_in), 0.5);
594 287064 : PetscScalar value_vec_tt = _TR * rho_old_interp * h_old_interp * volume / _dt;
595 287064 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
596 287064 : LibmeshPetscCall(
597 : VecSetValues(_hc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
598 :
599 : /// Advective derivative term
600 287064 : if (iz == first_node)
601 : {
602 12420 : PetscScalar value_vec_at = (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
603 12420 : PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
604 12420 : LibmeshPetscCall(VecSetValues(
605 : _hc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
606 : }
607 : else
608 : {
609 274644 : PetscInt row_at = i_ch + _n_channels * iz_ind;
610 274644 : PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
611 274644 : PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in);
612 274644 : LibmeshPetscCall(MatSetValues(
613 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
614 : }
615 :
616 : // Adding diagonal elements
617 287064 : PetscInt row_at = i_ch + _n_channels * iz_ind;
618 287064 : PetscInt col_at = i_ch + _n_channels * iz_ind;
619 287064 : PetscScalar value_at = (*_mdot_soln)(node_out);
620 287064 : LibmeshPetscCall(MatSetValues(
621 : _hc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
622 :
623 : /// Cross derivative term
624 : unsigned int counter = 0;
625 : unsigned int cross_index = iz; // iz-1;
626 1186488 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
627 : {
628 899424 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
629 : unsigned int ii_ch = chans.first;
630 : unsigned int jj_ch = chans.second;
631 899424 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
632 899424 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
633 : PetscScalar h_star;
634 : // figure out donor axial velocity
635 899424 : if (_Wij(i_gap, cross_index) > 0.0)
636 : {
637 529524 : if (iz == first_node)
638 : {
639 22320 : h_star = (*_h_soln)(node_in_i);
640 44640 : PetscScalar value_vec_ct = -1.0 * alpha *
641 22320 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
642 22320 : _Wij(i_gap, cross_index) * h_star;
643 22320 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
644 22320 : LibmeshPetscCall(VecSetValues(
645 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
646 : }
647 : else
648 : {
649 507204 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
650 507204 : _Wij(i_gap, cross_index);
651 507204 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
652 507204 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
653 507204 : LibmeshPetscCall(MatSetValues(
654 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
655 : }
656 529524 : PetscScalar value_ct = (1.0 - alpha) *
657 529524 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
658 529524 : _Wij(i_gap, cross_index);
659 529524 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
660 529524 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
661 529524 : LibmeshPetscCall(MatSetValues(
662 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
663 : }
664 369900 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
665 : {
666 319212 : if (iz == first_node)
667 : {
668 16128 : h_star = (*_h_soln)(node_in_j);
669 32256 : PetscScalar value_vec_ct = -1.0 * alpha *
670 16128 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
671 16128 : _Wij(i_gap, cross_index) * h_star;
672 16128 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
673 16128 : LibmeshPetscCall(VecSetValues(
674 : _hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
675 : }
676 : else
677 : {
678 303084 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
679 303084 : _Wij(i_gap, cross_index);
680 303084 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
681 303084 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
682 303084 : LibmeshPetscCall(MatSetValues(
683 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
684 : }
685 319212 : PetscScalar value_ct = (1.0 - alpha) *
686 319212 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
687 319212 : _Wij(i_gap, cross_index);
688 319212 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
689 319212 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
690 319212 : LibmeshPetscCall(MatSetValues(
691 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
692 : }
693 :
694 : // Turbulent cross flows
695 899424 : if (iz == first_node)
696 : {
697 : PetscScalar value_vec_ct =
698 39888 : -2.0 * alpha * (*_h_soln)(node_in)*_WijPrime(i_gap, cross_index);
699 39888 : value_vec_ct += alpha * (*_h_soln)(node_in_j)*_WijPrime(i_gap, cross_index);
700 39888 : value_vec_ct += alpha * (*_h_soln)(node_in_i)*_WijPrime(i_gap, cross_index);
701 39888 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
702 39888 : LibmeshPetscCall(
703 : VecSetValues(_hc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
704 : }
705 : else
706 : {
707 859536 : PetscScalar value_center_ct = 2.0 * alpha * _WijPrime(i_gap, cross_index);
708 859536 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
709 859536 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
710 859536 : LibmeshPetscCall(MatSetValues(
711 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
712 :
713 859536 : PetscScalar value_left_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
714 859536 : row_ct = i_ch + _n_channels * iz_ind;
715 859536 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
716 859536 : LibmeshPetscCall(MatSetValues(
717 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
718 :
719 859536 : PetscScalar value_right_ct = -1.0 * alpha * _WijPrime(i_gap, cross_index);
720 859536 : row_ct = i_ch + _n_channels * iz_ind;
721 859536 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
722 859536 : LibmeshPetscCall(MatSetValues(
723 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
724 : }
725 899424 : PetscScalar value_center_ct = 2.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
726 899424 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
727 899424 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
728 899424 : LibmeshPetscCall(MatSetValues(
729 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
730 :
731 899424 : PetscScalar value_left_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
732 899424 : row_ct = i_ch + _n_channels * iz_ind;
733 899424 : col_ct = jj_ch + _n_channels * iz_ind;
734 899424 : LibmeshPetscCall(MatSetValues(
735 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
736 :
737 899424 : PetscScalar value_right_ct = -1.0 * (1.0 - alpha) * _WijPrime(i_gap, cross_index);
738 899424 : row_ct = i_ch + _n_channels * iz_ind;
739 899424 : col_ct = ii_ch + _n_channels * iz_ind;
740 899424 : LibmeshPetscCall(MatSetValues(
741 : _hc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
742 899424 : counter++;
743 : }
744 :
745 : /// Added heat enthalpy
746 287064 : PetscScalar added_enthalpy = computeAddedHeatPin(i_ch, iz);
747 287064 : PetscInt row_vec_ht = i_ch + _n_channels * iz_ind;
748 287064 : LibmeshPetscCall(
749 : VecSetValues(_hc_added_heat_rhs, 1, &row_vec_ht, &added_enthalpy, ADD_VALUES));
750 : }
751 : }
752 : /// Assembling system
753 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
754 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
755 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
756 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
757 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
758 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
759 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
760 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
761 : // Matrix
762 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
763 474 : LibmeshPetscCall(MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
764 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
765 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
766 474 : LibmeshPetscCall(
767 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
768 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
769 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
770 474 : LibmeshPetscCall(
771 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
772 : #else
773 : LibmeshPetscCall(
774 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
775 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
776 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
777 : LibmeshPetscCall(
778 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
779 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
780 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
781 : LibmeshPetscCall(
782 : MatAXPY(_hc_sys_h_mat, 1.0, _hc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
783 : #endif
784 474 : LibmeshPetscCall(MatAssemblyBegin(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
785 474 : LibmeshPetscCall(MatAssemblyEnd(_hc_sys_h_mat, MAT_FINAL_ASSEMBLY));
786 474 : if (_verbose_subchannel)
787 294 : _console << "Block: " << iblock << " - Enthalpy conservation matrix assembled" << std::endl;
788 : // RHS
789 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_time_derivative_rhs));
790 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_advective_derivative_rhs));
791 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs));
792 474 : LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs));
793 :
794 474 : if (_segregated_bool || (!_monolithic_thermal_bool))
795 : {
796 : // Assembly the matrix system
797 : KSP ksploc;
798 : PC pc;
799 : Vec sol;
800 360 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
801 360 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
802 360 : LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
803 360 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
804 360 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
805 360 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
806 360 : LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_"));
807 360 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
808 360 : LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
809 : PetscScalar * xx;
810 360 : LibmeshPetscCall(VecGetArray(sol, &xx));
811 11496 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
812 : {
813 11136 : auto iz_ind = iz - first_node;
814 257160 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
815 : {
816 246024 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
817 246024 : auto h_out = xx[iz_ind * _n_channels + i_ch];
818 246024 : if (h_out < 0)
819 : {
820 0 : mooseError(name(),
821 : " : Calculation of negative Enthalpy h_out = : ",
822 : h_out,
823 : " Axial Level= : ",
824 : iz);
825 : }
826 246024 : _h_soln->set(node_out, h_out);
827 : }
828 : }
829 360 : LibmeshPetscCall(KSPDestroy(&ksploc));
830 360 : LibmeshPetscCall(VecDestroy(&sol));
831 : }
832 : }
833 668 : }
|
