12 #include "libmesh/petsc_vector.h" 13 #include "libmesh/dense_matrix.h" 14 #include "libmesh/dense_vector.h" 29 const PetscScalar * xx;
39 for (PetscInt i = 0; i < size; i++)
40 solution_seed(i) = xx[i];
48 for (
int i = 0; i < size; i++)
49 ff[i] = Wij_residual_vector(i);
58 MooseEnum schemes(
"upwind downwind central_difference exponential",
"central_difference");
62 params.
addRequiredParam<
unsigned int>(
"n_blocks",
"The number of blocks in the axial direction");
64 params.
addParam<
Real>(
"P_tol", 1e-6,
"Pressure tolerance");
65 params.
addParam<
Real>(
"T_tol", 1e-6,
"Temperature tolerance");
66 params.
addParam<
int>(
"T_maxit", 100,
"Maximum number of iterations for inner temperature loop");
67 params.
addParam<PetscReal>(
"rtol", 1e-6,
"Relative tolerance for ksp solver");
68 params.
addParam<PetscReal>(
"atol", 1e-6,
"Absolute tolerance for ksp solver");
69 params.
addParam<PetscReal>(
"dtol", 1e5,
"Divergence tolerance or ksp solver");
70 params.
addParam<PetscInt>(
"maxit", 1e4,
"Maximum number of iterations for ksp solver");
73 "Interpolation scheme used for the method. Default is exponential");
75 "implicit",
false,
"Boolean to define the use of explicit or implicit solution.");
77 "staggered_pressure",
false,
"Boolean to define the use of explicit or implicit solution.");
79 "segregated",
true,
"Boolean to define whether to use a segregated solution.");
81 "monolithic_thermal",
false,
"Boolean to define whether to use thermal monolithic solve.");
83 "verbose_subchannel",
false,
"Boolean to print out information related to subchannel solve.");
87 "Boolean that activates the deformation effect based on values for: displacement, Dpin");
88 params.
addRequiredParam<
bool>(
"compute_density",
"Flag that enables the calculation of density");
90 "Flag that enables the calculation of viscosity");
93 "Flag that informs whether we solve the Enthalpy/Temperature equations or not");
95 "P_out",
"The postprocessor (or scalar) that provides the value of outlet pressure [Pa]");
96 params.
addRequiredParam<UserObjectName>(
"fp",
"Fluid properties user object name");
104 _n_blocks(getParam<unsigned
int>(
"n_blocks")),
105 _Wij(declareRestartableData<
libMesh::DenseMatrix<
Real>>(
"Wij")),
107 _kij(_subchannel_mesh.getKij()),
109 _compute_density(getParam<bool>(
"compute_density")),
110 _compute_viscosity(getParam<bool>(
"compute_viscosity")),
111 _compute_power(getParam<bool>(
"compute_power")),
112 _pin_mesh_exist(_subchannel_mesh.pinMeshExist()),
113 _duct_mesh_exist(_subchannel_mesh.ductMeshExist()),
114 _P_out(getPostprocessorValue(
"P_out")),
115 _CT(getParam<
Real>(
"CT")),
116 _P_tol(getParam<
Real>(
"P_tol")),
117 _T_tol(getParam<
Real>(
"T_tol")),
118 _T_maxit(getParam<
int>(
"T_maxit")),
119 _rtol(getParam<PetscReal>(
"rtol")),
120 _atol(getParam<PetscReal>(
"atol")),
121 _dtol(getParam<PetscReal>(
"dtol")),
122 _maxit(getParam<PetscInt>(
"maxit")),
123 _interpolation_scheme(getParam<
MooseEnum>(
"interpolation_scheme")),
124 _implicit_bool(getParam<bool>(
"implicit")),
125 _staggered_pressure_bool(getParam<bool>(
"staggered_pressure")),
126 _segregated_bool(getParam<bool>(
"segregated")),
127 _monolithic_thermal_bool(getParam<bool>(
"monolithic_thermal")),
128 _verbose_subchannel(getParam<bool>(
"verbose_subchannel")),
129 _deformation(getParam<bool>(
"deformation")),
132 _q_prime_duct_soln(nullptr),
137 mooseError(
"Cannot use more than one MPI process.");
159 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
235 ": When implicit number of blocks can't be equal to number of cells. This will " 236 "cause problems with the subchannel interpolation scheme.");
244 _fp = &getUserObject<SinglePhaseFluidProperties>(getParam<UserObjectName>(
"fp"));
273 PetscErrorCode ierr =
cleanUp();
285 LibmeshPetscCall(VecDestroy(&
_Wij_vec));
286 LibmeshPetscCall(VecDestroy(&
_prod));
287 LibmeshPetscCall(VecDestroy(&
_prodp));
355 ": Interpolation scheme should be a string: upwind, downwind, central_difference, " 362 PetscScalar botValue,
366 return alpha * botValue + (1.0 -
alpha) * topValue;
372 unsigned int last_node = (iblock + 1) *
_block_size;
373 unsigned int first_node = iblock *
_block_size + 1;
376 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
378 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
380 int i =
_n_gaps * (iz - first_node) + i_gap;
381 solution_seed(i) =
_Wij(i_gap, iz);
392 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
394 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
405 unsigned int last_node = (iblock + 1) *
_block_size;
406 unsigned int first_node = iblock *
_block_size + 1;
410 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
412 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
417 unsigned int counter = 0;
431 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
433 unsigned int iz_ind = iz - first_node;
434 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
437 unsigned int counter = 0;
441 PetscInt col = i_gap +
_n_gaps * iz_ind;
448 LibmeshPetscCall(MatAssemblyBegin(
_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
449 LibmeshPetscCall(MatAssemblyEnd(
_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
455 LibmeshPetscCall(VecDuplicate(
_Wij_vec, &loc_Wij));
457 loc_Wij,
_Wij, first_node, last_node + 1,
_n_gaps));
459 LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
461 LibmeshPetscCall(VecDestroy(&loc_prod));
462 LibmeshPetscCall(VecDestroy(&loc_Wij));
470 unsigned int last_node = (iblock + 1) *
_block_size;
471 unsigned int first_node = iblock *
_block_size + 1;
474 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
477 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
481 auto volume = dz * (*_S_flow_soln)(node_in);
484 auto mdot_out = (*_mdot_soln)(node_in) - (*
_SumWij_soln)(node_out)-time_term;
489 " : Calculation of negative mass flow mdot_out = : ",
493 " - Implicit solves are required for recirculating flow.");
501 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
504 auto iz_ind = iz - first_node;
505 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
509 auto volume = dz * (*_S_flow_soln)(node_in);
514 PetscScalar value_vec = -1.0 * time_term;
519 if (iz == first_node)
521 PetscScalar value_vec = (*_mdot_soln)(node_in);
530 PetscScalar
value = -1.0;
538 PetscScalar
value = 1.0;
545 PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out);
561 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
563 LibmeshPetscCall(KSPGetPC(ksploc, &pc));
564 LibmeshPetscCall(PCSetType(pc, PCJACOBI));
566 LibmeshPetscCall(KSPSetFromOptions(ksploc));
568 LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
571 LibmeshPetscCall(KSPDestroy(&ksploc));
572 LibmeshPetscCall(VecDestroy(&sol));
580 unsigned int last_node = (iblock + 1) *
_block_size;
581 unsigned int first_node = iblock *
_block_size + 1;
584 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
588 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
592 auto rho_in = (*_rho_soln)(node_in);
593 auto rho_out = (*_rho_soln)(node_out);
594 auto mu_in = (*_mu_soln)(node_in);
595 auto S = (*_S_flow_soln)(node_in);
596 auto w_perim = (*_w_perim_soln)(node_in);
598 auto Dh_i = 4.0 *
S / w_perim;
599 auto time_term =
_TR * ((*_mdot_soln)(node_out)-
_mdot_soln->old(node_out)) * dz /
_dt -
604 auto mass_term2 = -2.0 * (*_mdot_soln)(node_out) * (*
_SumWij_soln)(node_out) /
S / rho_in;
605 auto crossflow_term = 0.0;
606 auto turbulent_term = 0.0;
607 unsigned int counter = 0;
611 unsigned int ii_ch = chans.first;
612 unsigned int jj_ch = chans.second;
617 auto rho_i = (*_rho_soln)(node_in_i);
618 auto rho_j = (*_rho_soln)(node_in_j);
619 auto Si = (*_S_flow_soln)(node_in_i);
620 auto Sj = (*_S_flow_soln)(node_in_j);
623 if (
_Wij(i_gap, iz) > 0.0)
624 u_star = (*_mdot_soln)(node_out_i) / Si / rho_i;
626 u_star = (*_mdot_soln)(node_out_j) / Sj / rho_j;
631 turbulent_term +=
_WijPrime(i_gap, iz) * (2 * (*_mdot_soln)(node_out) / rho_in /
S -
636 turbulent_term *=
_CT;
637 auto Re = (((*_mdot_soln)(node_in) /
S) * Dh_i / mu_in);
646 ki = k_grid[i_ch][iz - 1];
648 ki = k_grid[i_ch][iz];
649 auto friction_term = (fi * dz / Dh_i + ki) * 0.5 *
651 (
S * (*_rho_soln)(node_out));
652 auto gravity_term =
_g_grav * (*_rho_soln)(node_out)*dz *
S;
653 auto DP =
std::pow(
S, -1.0) * (time_term + mass_term1 + mass_term2 + crossflow_term +
654 turbulent_term + friction_term + gravity_term);
672 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
676 auto iz_ind = iz - first_node;
677 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
684 PetscScalar Pe = 0.5;
688 auto S_in = (*_S_flow_soln)(node_in);
689 auto S_out = (*_S_flow_soln)(node_out);
691 auto w_perim_in = (*_w_perim_soln)(node_in);
692 auto w_perim_out = (*_w_perim_soln)(node_out);
696 auto mu_in = (*_mu_soln)(node_in);
697 auto mu_out = (*_mu_soln)(node_out);
699 auto Dh_i = 4.0 * S_interp / w_perim_interp;
701 auto Re = ((mdot_loc / S_interp) * Dh_i / mu_interp);
710 ki = k_grid[i_ch][iz - 1];
712 ki = k_grid[i_ch][iz];
713 Pe = 1.0 / ((fi * dz / Dh_i + ki) * 0.5) * mdot_loc / std::abs(mdot_loc);
718 auto rho_in = (*_rho_soln)(node_in);
719 auto rho_out = (*_rho_soln)(node_out);
723 auto mu_in = (*_mu_soln)(node_in);
724 auto mu_out = (*_mu_soln)(node_out);
728 auto S_in = (*_S_flow_soln)(node_in);
729 auto S_out = (*_S_flow_soln)(node_out);
733 auto w_perim_in = (*_w_perim_soln)(node_in);
734 auto w_perim_out = (*_w_perim_soln)(node_out);
738 auto Dh_i = 4.0 * S_interp / w_perim_interp;
741 if (iz == first_node)
743 PetscScalar value_vec_tt = -1.0 *
_TR *
alpha * (*_mdot_soln)(node_in)*dz /
_dt;
751 PetscInt col_tt = i_ch +
_n_channels * (iz_ind - 1);
753 LibmeshPetscCall(MatSetValues(
760 PetscScalar value_tt =
_TR * (1.0 -
alpha) * dz /
_dt;
761 LibmeshPetscCall(MatSetValues(
765 PetscScalar mdot_old_interp =
767 PetscScalar value_vec_tt =
_TR * mdot_old_interp * dz /
_dt;
773 if (iz == first_node)
777 LibmeshPetscCall(VecSetValues(
783 PetscInt col_at = i_ch +
_n_channels * (iz_ind - 1);
784 PetscScalar value_at = -1.0 * std::abs((*
_mdot_soln)(node_in)) / (S_in * rho_in);
785 LibmeshPetscCall(MatSetValues(
792 PetscScalar value_at = std::abs((*
_mdot_soln)(node_out)) / (S_out * rho_out);
793 LibmeshPetscCall(MatSetValues(
797 unsigned int counter = 0;
798 unsigned int cross_index = iz;
802 unsigned int ii_ch = chans.first;
803 unsigned int jj_ch = chans.second;
818 if (
_Wij(i_gap, cross_index) > 0.0)
820 if (iz == first_node)
822 u_star = (*_mdot_soln)(node_in_i) / S_i / rho_i;
823 PetscScalar value_vec_ct = -1.0 *
alpha *
825 _Wij(i_gap, cross_index) * u_star;
827 LibmeshPetscCall(VecSetValues(
833 _Wij(i_gap, cross_index) / S_i / rho_i;
835 PetscInt col_ct = ii_ch +
_n_channels * (iz_ind - 1);
836 LibmeshPetscCall(MatSetValues(
839 PetscScalar value_ct = (1.0 -
alpha) *
841 _Wij(i_gap, cross_index) / S_i / rho_i;
844 LibmeshPetscCall(MatSetValues(
847 else if (
_Wij(i_gap, cross_index) < 0.0)
849 if (iz == first_node)
851 u_star = (*_mdot_soln)(node_in_j) / S_j / rho_j;
852 PetscScalar value_vec_ct = -1.0 *
alpha *
854 _Wij(i_gap, cross_index) * u_star;
856 LibmeshPetscCall(VecSetValues(
862 _Wij(i_gap, cross_index) / S_j / rho_j;
864 PetscInt col_ct = jj_ch +
_n_channels * (iz_ind - 1);
865 LibmeshPetscCall(MatSetValues(
868 PetscScalar value_ct = (1.0 -
alpha) *
870 _Wij(i_gap, cross_index) / S_j / rho_j;
873 LibmeshPetscCall(MatSetValues(
877 if (iz == first_node)
879 PetscScalar value_vec_ct = -2.0 *
alpha * (*_mdot_soln)(node_in)*
_CT *
880 _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
881 value_vec_ct +=
alpha * (*_mdot_soln)(node_in_j)*
_CT *
_WijPrime(i_gap, cross_index) /
883 value_vec_ct +=
alpha * (*_mdot_soln)(node_in_i)*
_CT *
_WijPrime(i_gap, cross_index) /
891 PetscScalar value_center_ct =
894 PetscInt col_ct = i_ch +
_n_channels * (iz_ind - 1);
895 LibmeshPetscCall(MatSetValues(
898 PetscScalar value_left_ct =
902 LibmeshPetscCall(MatSetValues(
905 PetscScalar value_right_ct =
909 LibmeshPetscCall(MatSetValues(
913 PetscScalar value_center_ct =
914 2.0 * (1.0 -
alpha) *
_CT *
_WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
917 LibmeshPetscCall(MatSetValues(
920 PetscScalar value_left_ct =
924 LibmeshPetscCall(MatSetValues(
927 PetscScalar value_right_ct =
931 LibmeshPetscCall(MatSetValues(
937 PetscScalar mdot_interp =
939 auto Re = ((mdot_interp / S_interp) * Dh_i / mu_interp);
948 ki = k_grid[i_ch][iz - 1];
950 ki = k_grid[i_ch][iz];
951 auto coef = (fi * dz / Dh_i + ki) * 0.5 * std::abs((*
_mdot_soln)(node_out)) /
952 (S_interp * rho_interp);
953 if (iz == first_node)
955 PetscScalar value_vec = -1.0 *
alpha * coef * (*_mdot_soln)(node_in);
977 PetscScalar value_vec = -1.0 *
_g_grav * rho_interp * dz * S_interp;
979 LibmeshPetscCall(VecSetValues(
_amc_gravity_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
996 #if !PETSC_VERSION_LESS_THAN(3, 15, 0) 1038 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1045 LibmeshPetscCall(VecGetArray(ls, &xx));
1046 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1048 auto iz_ind = iz - first_node;
1049 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1056 auto S_in = (*_S_flow_soln)(node_in);
1057 auto S_out = (*_S_flow_soln)(node_out);
1073 LibmeshPetscCall(VecDestroy(&ls));
1081 unsigned int last_node = (iblock + 1) *
_block_size;
1082 unsigned int first_node = iblock *
_block_size + 1;
1087 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1090 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1101 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1104 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1112 PetscScalar Pe = 0.5;
1114 if (iz == last_node)
1133 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1135 auto iz_ind = iz - first_node;
1137 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1143 auto S_in = (*_S_flow_soln)(node_in);
1144 auto S_out = (*_S_flow_soln)(node_out);
1150 PetscScalar
value = -1.0 * S_interp;
1154 if (iz == last_node)
1156 PetscScalar
value = -1.0 * (*_P_soln)(node_out)*S_interp;
1164 PetscScalar
value = 1.0 * S_interp;
1171 auto dp_out = (*_DP_soln)(node_out);
1172 PetscScalar value_v = -1.0 * dp_out * S_interp;
1188 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1190 LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1191 LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1193 LibmeshPetscCall(KSPSetFromOptions(ksploc));
1196 LibmeshPetscCall(VecGetArray(sol, &xx));
1198 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1200 auto iz_ind = iz - first_node;
1201 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1209 LibmeshPetscCall(KSPDestroy(&ksploc));
1210 LibmeshPetscCall(VecDestroy(&sol));
1216 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1218 auto iz_ind = iz - first_node;
1220 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1226 auto S_in = (*_S_flow_soln)(node_in);
1227 auto S_out = (*_S_flow_soln)(node_out);
1233 PetscScalar
value = -1.0 * S_interp;
1237 if (iz == last_node)
1239 PetscScalar
value = -1.0 * (*_P_soln)(node_out)*S_interp;
1243 auto dp_out = (*_DP_soln)(node_out);
1244 PetscScalar value_v = -1.0 * dp_out / 2.0 * S_interp;
1253 PetscScalar
value = 1.0 * S_interp;
1259 auto dp_in = (*_DP_soln)(node_in);
1260 auto dp_out = (*_DP_soln)(node_out);
1262 PetscScalar value_v = -1.0 * dp_interp * S_interp;
1274 _console <<
"Block: " << iblock <<
" - Axial momentum pressure force matrix assembled" 1283 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1285 LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1286 LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1288 LibmeshPetscCall(KSPSetFromOptions(ksploc));
1291 LibmeshPetscCall(VecGetArray(sol, &xx));
1293 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
1295 auto iz_ind = iz - first_node;
1296 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1304 LibmeshPetscCall(KSPDestroy(&ksploc));
1305 LibmeshPetscCall(VecDestroy(&sol));
1324 auto heat_rate_in = (*_q_prime_duct_soln)(node_in_duct);
1325 auto heat_rate_out = (*_q_prime_duct_soln)(node_out_duct);
1342 unsigned int last_node = (iblock + 1) *
_block_size;
1343 unsigned int first_node = iblock *
_block_size + 1;
1344 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1346 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1357 unsigned int last_node = (iblock + 1) *
_block_size;
1358 unsigned int first_node = iblock *
_block_size + 1;
1361 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1367 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1369 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1380 unsigned int last_node = (iblock + 1) *
_block_size;
1381 unsigned int first_node = iblock *
_block_size + 1;
1384 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1390 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1392 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1406 unsigned int last_node = (iblock + 1) *
_block_size;
1409 for (
unsigned int iz = first_node + 1; iz < last_node + 1; iz++)
1412 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1415 unsigned int i_ch = chans.first;
1416 unsigned int j_ch = chans.second;
1421 auto rho_i = (*_rho_soln)(node_in_i);
1422 auto rho_j = (*_rho_soln)(node_in_j);
1423 auto Si = (*_S_flow_soln)(node_in_i);
1424 auto Sj = (*_S_flow_soln)(node_in_j);
1428 auto friction_term =
_kij *
_Wij(i_gap, iz) * std::abs(
_Wij(i_gap, iz));
1429 auto DPij = (*_P_soln)(node_in_i) - (*
_P_soln)(node_in_j);
1431 auto rho_star = 0.0;
1432 if (
_Wij(i_gap, iz) > 0.0)
1434 else if (
_Wij(i_gap, iz) < 0.0)
1437 rho_star = (rho_i + rho_j) / 2.0;
1438 auto mass_term_out =
1439 (*_mdot_soln)(node_out_i) / Si / rho_i + (*
_mdot_soln)(node_out_j) / Sj / rho_j;
1441 (*_mdot_soln)(node_in_i) / Si / rho_i + (*
_mdot_soln)(node_in_j) / Sj / rho_j;
1442 auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out *
_Wij(i_gap, iz);
1443 auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in *
_Wij(i_gap, iz - 1);
1444 auto inertia_term = term_out - term_in;
1445 auto pressure_term = 2 *
std::pow(Sij, 2.0) * DPij * rho_star;
1450 time_term + friction_term + inertia_term - pressure_term;
1467 unsigned int last_node = (iblock + 1) *
_block_size;
1470 for (
unsigned int iz = first_node + 1; iz < last_node + 1; iz++)
1473 auto iz_ind = iz - first_node - 1;
1474 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1477 unsigned int i_ch = chans.first;
1478 unsigned int j_ch = chans.second;
1485 auto rho_i_in = (*_rho_soln)(node_in_i);
1486 auto rho_i_out = (*_rho_soln)(node_out_i);
1488 auto rho_j_in = (*_rho_soln)(node_in_j);
1489 auto rho_j_out = (*_rho_soln)(node_out_j);
1493 auto S_i_in = (*_S_flow_soln)(node_in_i);
1494 auto S_i_out = (*_S_flow_soln)(node_out_i);
1495 auto S_j_in = (*_S_flow_soln)(node_in_j);
1496 auto S_j_out = (*_S_flow_soln)(node_out_j);
1503 auto rho_star = 0.0;
1504 if (
_Wij(i_gap, iz) > 0.0)
1505 rho_star = rho_i_interp;
1506 else if (
_Wij(i_gap, iz) < 0.0)
1507 rho_star = rho_j_interp;
1509 rho_star = (rho_i_interp + rho_j_interp) / 2.0;
1512 PetscScalar time_factor =
_TR * Lij * Sij * rho_star /
_dt;
1513 PetscInt row_td = i_gap +
_n_gaps * iz_ind;
1514 PetscInt col_td = i_gap +
_n_gaps * iz_ind;
1515 PetscScalar value_td = time_factor;
1516 LibmeshPetscCall(MatSetValues(
1518 PetscScalar value_td_rhs = time_factor *
_Wij_old(i_gap, iz);
1523 PetscScalar Pe = 0.5;
1525 auto mass_term_out = (*_mdot_soln)(node_out_i) / S_i_out / rho_i_out +
1526 (*
_mdot_soln)(node_out_j) / S_j_out / rho_j_out;
1527 auto mass_term_in = (*_mdot_soln)(node_in_i) / S_i_in / rho_i_in +
1528 (*
_mdot_soln)(node_in_j) / S_j_in / rho_j_in;
1529 auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out / 2.0;
1530 auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in / 2.0;
1531 if (iz == first_node + 1)
1533 PetscInt row_ad = i_gap +
_n_gaps * iz_ind;
1534 PetscScalar value_ad = term_in *
alpha *
_Wij(i_gap, iz - 1);
1538 PetscInt col_ad = i_gap +
_n_gaps * iz_ind;
1539 value_ad = -1.0 * term_in * (1.0 -
alpha) + term_out *
alpha;
1540 LibmeshPetscCall(MatSetValues(
1543 col_ad = i_gap +
_n_gaps * (iz_ind + 1);
1544 value_ad = term_out * (1.0 -
alpha);
1545 LibmeshPetscCall(MatSetValues(
1548 else if (iz == last_node)
1550 PetscInt row_ad = i_gap +
_n_gaps * iz_ind;
1551 PetscInt col_ad = i_gap +
_n_gaps * (iz_ind - 1);
1552 PetscScalar value_ad = -1.0 * term_in *
alpha;
1553 LibmeshPetscCall(MatSetValues(
1556 col_ad = i_gap +
_n_gaps * iz_ind;
1557 value_ad = -1.0 * term_in * (1.0 -
alpha) + term_out *
alpha;
1558 LibmeshPetscCall(MatSetValues(
1561 value_ad = -1.0 * term_out * (1.0 -
alpha) *
_Wij(i_gap, iz);
1567 PetscInt row_ad = i_gap +
_n_gaps * iz_ind;
1568 PetscInt col_ad = i_gap +
_n_gaps * (iz_ind - 1);
1569 PetscScalar value_ad = -1.0 * term_in *
alpha;
1570 LibmeshPetscCall(MatSetValues(
1573 col_ad = i_gap +
_n_gaps * iz_ind;
1574 value_ad = -1.0 * term_in * (1.0 -
alpha) + term_out *
alpha;
1575 LibmeshPetscCall(MatSetValues(
1578 col_ad = i_gap +
_n_gaps * (iz_ind + 1);
1579 value_ad = term_out * (1.0 -
alpha);
1580 LibmeshPetscCall(MatSetValues(
1584 PetscInt row_ff = i_gap +
_n_gaps * iz_ind;
1585 PetscInt col_ff = i_gap +
_n_gaps * iz_ind;
1586 PetscScalar value_ff =
_kij * std::abs(
_Wij(i_gap, iz)) / 2.0;
1587 LibmeshPetscCall(MatSetValues(
1595 PetscScalar pressure_factor =
std::pow(Sij, 2.0) * rho_star;
1596 PetscInt row_pf = i_gap +
_n_gaps * iz_ind;
1598 PetscScalar value_pf = -1.0 *
alpha * pressure_factor;
1602 value_pf =
alpha * pressure_factor;
1606 if (iz == last_node)
1608 PetscInt row_pf = i_gap +
_n_gaps * iz_ind;
1609 PetscScalar value_pf = (1.0 -
alpha) * pressure_factor * (*
_P_soln)(node_out_i);
1612 value_pf = -1.0 * (1.0 -
alpha) * pressure_factor * (*
_P_soln)(node_out_j);
1618 row_pf = i_gap +
_n_gaps * iz_ind;
1620 value_pf = -1.0 * (1.0 -
alpha) * pressure_factor;
1621 LibmeshPetscCall(MatSetValues(
1624 value_pf = (1.0 -
alpha) * pressure_factor;
1625 LibmeshPetscCall(MatSetValues(
1631 PetscScalar pressure_factor =
std::pow(Sij, 2.0) * rho_star;
1632 PetscInt row_pf = i_gap +
_n_gaps * iz_ind;
1634 PetscScalar value_pf = -1.0 * pressure_factor;
1638 value_pf = pressure_factor;
1658 #if !PETSC_VERSION_LESS_THAN(3, 15, 0) 1695 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1701 LibmeshPetscCall(VecAXPY(sol_holder_W, 1.0, sol_holder_P));
1703 LibmeshPetscCall(VecGetArray(sol_holder_W, &xx));
1704 for (
unsigned int iz = first_node + 1; iz < last_node + 1; iz++)
1706 auto iz_ind = iz - first_node - 1;
1707 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1712 LibmeshPetscCall(VecDestroy(&sol_holder_P));
1713 LibmeshPetscCall(VecDestroy(&sol_holder_W));
1721 unsigned int last_node = (iblock + 1) *
_block_size;
1722 unsigned int first_node = iblock *
_block_size + 1;
1723 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1726 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1729 unsigned int i_ch = chans.first;
1730 unsigned int j_ch = chans.second;
1735 auto Si_in = (*_S_flow_soln)(node_in_i);
1736 auto Sj_in = (*_S_flow_soln)(node_in_j);
1737 auto Si_out = (*_S_flow_soln)(node_out_i);
1738 auto Sj_out = (*_S_flow_soln)(node_out_j);
1740 auto Sij = dz * gap;
1742 0.5 * (((*_mdot_soln)(node_in_i) + (*
_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
1748 _WijPrime(i_gap, iz) = beta * avg_massflux * Sij;
1752 auto iz_ind = iz - first_node;
1753 PetscScalar base_value = beta * 0.5 * Sij;
1756 if (iz == first_node)
1758 PetscScalar value_tl = -1.0 * base_value / (Si_in + Sj_in) *
1760 PetscInt row = i_gap +
_n_gaps * iz_ind;
1766 PetscScalar value_tl = base_value / (Si_in + Sj_in);
1767 PetscInt row = i_gap +
_n_gaps * iz_ind;
1769 PetscInt col_ich = i_ch +
_n_channels * (iz_ind - 1);
1770 LibmeshPetscCall(MatSetValues(
1773 PetscInt col_jch = j_ch +
_n_channels * (iz_ind - 1);
1774 LibmeshPetscCall(MatSetValues(
1779 PetscScalar value_bl = base_value / (Si_out + Sj_out);
1780 PetscInt row = i_gap +
_n_gaps * iz_ind;
1783 LibmeshPetscCall(MatSetValues(
1787 LibmeshPetscCall(MatSetValues(
1802 LibmeshPetscCall(VecDuplicate(
_Wij_vec, &loc_Wij));
1803 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1809 LibmeshPetscCall(VecDestroy(&loc_prod));
1810 LibmeshPetscCall(VecDestroy(&loc_Wij));
1817 unsigned int last_node = (iblock + 1) *
_block_size;
1822 for (
unsigned int iz = first_node + 1; iz < last_node + 1; iz++)
1824 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1826 _Wij(i_gap, iz) = solution(i);
1847 for (
unsigned int i_gap = 0; i_gap <
_n_gaps; i_gap++)
1853 return Wij_residual_vector;
1869 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1871 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP,
"Example is only for sequential runs");
1872 LibmeshPetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
1873 LibmeshPetscCall(VecCreate(PETSC_COMM_WORLD, &
x));
1875 LibmeshPetscCall(VecSetFromOptions(
x));
1876 LibmeshPetscCall(VecDuplicate(
x, &r));
1878 #if PETSC_VERSION_LESS_THAN(3, 13, 0) 1879 LibmeshPetscCall(PetscOptionsSetValue(PETSC_NULL,
"-snes_mf", PETSC_NULL));
1881 LibmeshPetscCall(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
1884 ctx.iblock = iblock;
1887 LibmeshPetscCall(SNESGetKSP(snes, &ksp));
1888 LibmeshPetscCall(KSPGetPC(ksp, &pc));
1889 LibmeshPetscCall(PCSetType(pc, PCNONE));
1891 LibmeshPetscCall(SNESSetFromOptions(snes));
1892 LibmeshPetscCall(VecGetArray(
x, &xx));
1895 xx[i] = solution(i);
1897 LibmeshPetscCall(VecRestoreArray(
x, &xx));
1899 LibmeshPetscCall(SNESSolve(snes, NULL,
x));
1900 LibmeshPetscCall(VecGetArray(
x, &xx));
1904 LibmeshPetscCall(VecRestoreArray(
x, &xx));
1905 LibmeshPetscCall(VecDestroy(&
x));
1906 LibmeshPetscCall(VecDestroy(&r));
1907 LibmeshPetscCall(SNESDestroy(&snes));
1914 bool lag_block_thermal_solve =
true;
1922 std::vector<Mat> mat_array(Q * Q);
1923 std::vector<Vec> vec_array(Q);
1926 bool _axial_mass_flow_tight_coupling =
true;
1927 bool _pressure_axial_momentum_tight_coupling =
true;
1928 bool _pressure_cross_momentum_tight_coupling =
true;
1929 unsigned int first_node = iblock *
_block_size + 1;
1930 unsigned int last_node = (iblock + 1) *
_block_size;
1951 _console <<
"Starting nested system." << std::endl;
1952 _console <<
"Number of simultaneous variables: " << Q << std::endl;
1955 PetscInt field_num = 0;
1958 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
1959 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
1960 mat_array[Q * field_num + 1] = NULL;
1961 if (_axial_mass_flow_tight_coupling)
1963 LibmeshPetscCall(MatDuplicate(
_mc_sumWij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2]));
1964 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
1965 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
1969 mat_array[Q * field_num + 2] = NULL;
1973 mat_array[Q * field_num + 3] = NULL;
1977 if (!_axial_mass_flow_tight_coupling)
1981 LibmeshPetscCall(VecSet(sumWij_loc, 0.0));
1982 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
1984 auto iz_ind = iz - first_node;
1985 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
1989 PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out);
1991 LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec_2, &value_vec_2, ADD_VALUES));
1994 LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sumWij_loc));
1995 LibmeshPetscCall(VecDestroy(&sumWij_loc));
1999 _console <<
"Mass ok." << std::endl;
2002 if (_pressure_axial_momentum_tight_coupling)
2005 MatDuplicate(
_amc_sys_mdot_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 0]));
2006 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
2007 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
2011 mat_array[Q * field_num + 0] = NULL;
2015 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2016 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2017 mat_array[Q * field_num + 2] = NULL;
2020 mat_array[Q * field_num + 3] = NULL;
2024 if (_pressure_axial_momentum_tight_coupling)
2030 unsigned int last_node = (iblock + 1) *
_block_size;
2031 unsigned int first_node = iblock *
_block_size + 1;
2032 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2038 LibmeshPetscCall(VecAXPY(vec_array[field_num], -1.0, ls));
2039 LibmeshPetscCall(VecDestroy(&ls));
2043 _console <<
"Lin mom OK." << std::endl;
2047 mat_array[Q * field_num + 0] = NULL;
2048 if (_pressure_cross_momentum_tight_coupling)
2052 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2053 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2057 mat_array[Q * field_num + 1] = NULL;
2060 LibmeshPetscCall(MatDuplicate(
_cmc_sys_Wij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2]));
2061 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
2062 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
2065 mat_array[Q * field_num + 3] = NULL;
2070 if (_pressure_cross_momentum_tight_coupling)
2078 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2083 LibmeshPetscCall(VecScale(sol_holder_P, 1.0));
2084 LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sol_holder_P));
2085 LibmeshPetscCall(VecDestroy(&sol_holder_P));
2089 _console <<
"Cross mom ok." << std::endl;
2095 mat_array[Q * field_num + 0] = NULL;
2096 mat_array[Q * field_num + 1] = NULL;
2097 mat_array[Q * field_num + 2] = NULL;
2098 LibmeshPetscCall(MatDuplicate(
_hc_sys_h_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 3]));
2099 if (lag_block_thermal_solve)
2101 LibmeshPetscCall(MatZeroEntries(mat_array[Q * field_num + 3]));
2102 LibmeshPetscCall(MatShift(mat_array[Q * field_num + 3], 1.0));
2104 LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
2105 LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
2106 LibmeshPetscCall(VecDuplicate(
_hc_sys_h_rhs, &vec_array[field_num]));
2107 LibmeshPetscCall(VecCopy(
_hc_sys_h_rhs, vec_array[field_num]));
2108 if (lag_block_thermal_solve)
2110 LibmeshPetscCall(VecZeroEntries(vec_array[field_num]));
2111 LibmeshPetscCall(VecShift(vec_array[field_num], 1.0));
2115 _console <<
"Energy ok." << std::endl;
2123 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2133 LibmeshPetscCall(VecSet(unity_vec, 1.0));
2142 LibmeshPetscCall(VecSet(unity_vec_Wij, 1.0));
2145 Vec _Wij_old_loc_vec;
2147 LibmeshPetscCall(MatMult(mat_array[Q],
_prod, mdot_estimate));
2148 LibmeshPetscCall(MatGetDiagonal(mat_array[Q + 1], pmat_diag));
2149 LibmeshPetscCall(VecAXPY(pmat_diag, 1e-10, unity_vec));
2150 LibmeshPetscCall(VecPointwiseDivide(p_estimate, mdot_estimate, pmat_diag));
2151 LibmeshPetscCall(MatMult(mat_array[2 * Q + 1], p_estimate, sol_holder_P));
2153 LibmeshPetscCall(MatGetDiagonal(mat_array[2 * Q + 2], diag_Wij_loc));
2154 LibmeshPetscCall(VecAXPY(diag_Wij_loc, 1e-10, unity_vec_Wij));
2155 LibmeshPetscCall(VecPointwiseDivide(Wij_estimate, sol_holder_P, diag_Wij_loc));
2158 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
2160 auto iz_ind = iz - first_node;
2161 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
2163 PetscScalar sumWij = 0.0;
2164 unsigned int counter = 0;
2168 unsigned int i_ch_loc = chans.first;
2169 PetscInt row_vec = i_ch_loc +
_n_channels * iz_ind;
2170 PetscScalar loc_Wij_value;
2171 LibmeshPetscCall(VecGetValues(sol_holder_P, 1, &row_vec, &loc_Wij_value));
2176 LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec, &sumWij, INSERT_VALUES));
2180 PetscScalar min_mdot;
2181 LibmeshPetscCall(VecAbs(
_prod));
2182 LibmeshPetscCall(VecMin(
_prod, NULL, &min_mdot));
2184 _console <<
"Minimum estimated mdot: " << min_mdot << std::endl;
2186 LibmeshPetscCall(VecAbs(sumWij_loc));
2187 LibmeshPetscCall(VecMax(sumWij_loc, NULL, &
_max_sumWij));
2193 _Wij_loc_vec,
_Wij, first_node, last_node,
_n_gaps));
2194 LibmeshPetscCall(VecAbs(_Wij_loc_vec));
2197 LibmeshPetscCall(VecAbs(_Wij_old_loc_vec));
2198 LibmeshPetscCall(VecAXPY(_Wij_loc_vec, -1.0, _Wij_old_loc_vec));
2199 PetscScalar relax_factor;
2200 LibmeshPetscCall(VecAbs(_Wij_loc_vec));
2201 #if !PETSC_VERSION_LESS_THAN(3, 16, 0) 2202 LibmeshPetscCall(VecMean(_Wij_loc_vec, &relax_factor));
2204 VecSum(_Wij_loc_vec, &relax_factor);
2209 _console <<
"Relax base value: " << relax_factor << std::endl;
2211 PetscScalar resistance_relaxation = 0.9;
2219 if (_added_K < 10 && _added_K >= 1.0)
2221 if (_added_K < 1.0 && _added_K >= 0.1)
2223 if (_added_K < 0.1 && _added_K >= 0.01)
2225 if (_added_K < 1e-2 && _added_K >= 1e-3)
2231 LibmeshPetscCall(VecScale(unity_vec_Wij,
_added_K));
2235 LibmeshPetscCall(MatDiagonalSet(mat_array[2 * Q + 2], unity_vec_Wij, ADD_VALUES));
2236 LibmeshPetscCall(VecDestroy(&mdot_estimate));
2237 LibmeshPetscCall(VecDestroy(&pmat_diag));
2238 LibmeshPetscCall(VecDestroy(&unity_vec));
2239 LibmeshPetscCall(VecDestroy(&p_estimate));
2240 LibmeshPetscCall(VecDestroy(&sol_holder_P));
2241 LibmeshPetscCall(VecDestroy(&diag_Wij_loc));
2242 LibmeshPetscCall(VecDestroy(&unity_vec_Wij));
2243 LibmeshPetscCall(VecDestroy(&Wij_estimate));
2244 LibmeshPetscCall(VecDestroy(&sumWij_loc));
2245 LibmeshPetscCall(VecDestroy(&_Wij_loc_vec));
2246 LibmeshPetscCall(VecDestroy(&_Wij_old_loc_vec));
2249 PetscScalar relaxation_factor_mdot, relaxation_factor_P, relaxation_factor_Wij;
2250 relaxation_factor_mdot = 1.0;
2251 relaxation_factor_P = 1.0;
2252 relaxation_factor_Wij = 0.1;
2256 _console <<
"Relax mdot: " << relaxation_factor_mdot << std::endl;
2257 _console <<
"Relax P: " << relaxation_factor_P << std::endl;
2258 _console <<
"Relax Wij: " << relaxation_factor_Wij << std::endl;
2261 PetscInt field_num = 0;
2263 LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_mdot));
2264 LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_mdot));
2265 LibmeshPetscCall(VecScale(diag_mdot, 1.0 / relaxation_factor_mdot));
2267 MatDiagonalSet(mat_array[Q * field_num + field_num], diag_mdot, INSERT_VALUES));
2268 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2270 LibmeshPetscCall(VecScale(diag_mdot, (1.0 - relaxation_factor_mdot)));
2271 LibmeshPetscCall(VecPointwiseMult(
_prod,
_prod, diag_mdot));
2272 LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0,
_prod));
2273 LibmeshPetscCall(VecDestroy(&diag_mdot));
2276 _console <<
"mdot relaxed" << std::endl;
2280 LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_P));
2281 LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_P));
2282 LibmeshPetscCall(VecScale(diag_P, 1.0 / relaxation_factor_P));
2283 LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_P, INSERT_VALUES));
2285 _console <<
"Mat assembled" << std::endl;
2286 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2288 LibmeshPetscCall(VecScale(diag_P, (1.0 - relaxation_factor_P)));
2289 LibmeshPetscCall(VecPointwiseMult(
_prod,
_prod, diag_P));
2290 LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0,
_prod));
2291 LibmeshPetscCall(VecDestroy(&diag_P));
2294 _console <<
"P relaxed" << std::endl;
2298 LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_Wij));
2299 LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_Wij));
2300 LibmeshPetscCall(VecScale(diag_Wij, 1.0 / relaxation_factor_Wij));
2301 LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_Wij, INSERT_VALUES));
2304 LibmeshPetscCall(VecScale(diag_Wij, (1.0 - relaxation_factor_Wij)));
2306 LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0,
_Wij_vec));
2307 LibmeshPetscCall(VecDestroy(&diag_Wij));
2310 _console <<
"Wij relaxed" << std::endl;
2313 _console <<
"Linear solver relaxed." << std::endl;
2316 LibmeshPetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, mat_array.data(), &A_nest));
2317 LibmeshPetscCall(VecCreateNest(PETSC_COMM_WORLD, Q, NULL, vec_array.data(), &b_nest));
2319 _console <<
"Nested system created." << std::endl;
2323 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
2324 LibmeshPetscCall(KSPSetType(ksp, KSPFGMRES));
2326 LibmeshPetscCall(KSPSetOperators(ksp, A_nest, A_nest));
2328 LibmeshPetscCall(KSPGetPC(ksp, &pc));
2329 LibmeshPetscCall(PCSetType(pc, PCFIELDSPLIT));
2332 std::vector<IS> rows(Q);
2335 LibmeshPetscCall(MatNestGetISs(A_nest, rows.data(), NULL));
2336 for (PetscInt
j = 0;
j < Q; ++
j)
2339 LibmeshPetscCall(ISDuplicate(rows[M], &expand1));
2341 LibmeshPetscCall(PCFieldSplitSetIS(pc, NULL, expand1));
2342 LibmeshPetscCall(ISDestroy(&expand1));
2345 _console <<
"Linear solver assembled." << std::endl;
2348 LibmeshPetscCall(VecDuplicate(b_nest, &x_nest));
2349 LibmeshPetscCall(VecSet(x_nest, 0.0));
2350 LibmeshPetscCall(KSPSolve(ksp, b_nest, x_nest));
2353 LibmeshPetscCall(VecDestroy(&b_nest));
2354 LibmeshPetscCall(MatDestroy(&A_nest));
2355 LibmeshPetscCall(KSPDestroy(&ksp));
2356 for (PetscInt i = 0; i < Q * Q; i++)
2358 LibmeshPetscCall(MatDestroy(&mat_array[i]));
2360 for (PetscInt i = 0; i < Q; i++)
2362 LibmeshPetscCall(VecDestroy(&vec_array[i]));
2365 _console <<
"Solver elements destroyed." << std::endl;
2368 Vec sol_mdot, sol_p, sol_Wij;
2370 _console <<
"Vectors created." << std::endl;
2373 LibmeshPetscCall(VecNestGetSubVecs(x_nest, &num_vecs, &loc_vecs));
2375 _console <<
"Starting extraction." << std::endl;
2377 LibmeshPetscCall(VecCopy(loc_vecs[0], sol_mdot));
2379 LibmeshPetscCall(VecCopy(loc_vecs[1], sol_p));
2381 LibmeshPetscCall(VecCopy(loc_vecs[2], sol_Wij));
2383 _console <<
"Getting individual field solutions from coupled solver." << std::endl;
2386 PetscScalar * sol_p_array;
2387 LibmeshPetscCall(VecGetArray(sol_p, &sol_p_array));
2388 PetscScalar * sol_Wij_array;
2389 LibmeshPetscCall(VecGetArray(sol_Wij, &sol_Wij_array));
2392 LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
2396 for (
unsigned int iz = last_node; iz > first_node - 1; iz--)
2398 auto iz_ind = iz - first_node;
2399 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
2409 sol_Wij,
_Wij, first_node, last_node - 1,
_n_gaps));
2414 if (lag_block_thermal_solve)
2420 LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
2422 LibmeshPetscCall(KSPGetPC(ksploc, &pc));
2423 LibmeshPetscCall(PCSetType(pc, PCJACOBI));
2425 LibmeshPetscCall(KSPSetFromOptions(ksploc));
2428 LibmeshPetscCall(VecGetArray(sol, &xx));
2429 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
2431 auto iz_ind = iz - first_node;
2432 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
2439 " : Calculation of negative Enthalpy h_out = : ",
2444 _h_soln->set(node_out, h_out);
2447 LibmeshPetscCall(KSPDestroy(&ksploc));
2448 LibmeshPetscCall(VecDestroy(&sol));
2454 LibmeshPetscCall(VecCopy(loc_vecs[3], sol_h));
2455 PetscScalar * sol_h_array;
2456 LibmeshPetscCall(VecGetArray(sol_h, &sol_h_array));
2457 for (
unsigned int iz = first_node; iz < last_node + 1; iz++)
2459 auto iz_ind = iz - first_node;
2460 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
2463 auto h_out = sol_h_array[iz_ind *
_n_channels + i_ch];
2467 " : Calculation of negative Enthalpy h_out = : ",
2472 _h_soln->set(node_out, h_out);
2475 LibmeshPetscCall(VecDestroy(&sol_h));
2481 LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
2486 LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2489 LibmeshPetscCall(VecAbs(
_prod));
2497 _console <<
"Solutions assigned to MOOSE variables." << std::endl;
2500 LibmeshPetscCall(VecDestroy(&x_nest));
2501 LibmeshPetscCall(VecDestroy(&sol_mdot));
2502 LibmeshPetscCall(VecDestroy(&sol_p));
2503 LibmeshPetscCall(VecDestroy(&sol_Wij));
2504 LibmeshPetscCall(VecDestroy(&sumWij_loc));
2506 _console <<
"Solutions destroyed." << std::endl;
2514 _console <<
"Executing subchannel solver\n";
2519 _console <<
"Solution initialized" << std::endl;
2521 unsigned int P_it = 0;
2522 unsigned int P_it_max;
2532 while ((P_error >
_P_tol && P_it < P_it_max))
2537 _console <<
"Reached maximum number of axial pressure iterations" << std::endl;
2540 _console <<
"Solving Outer Iteration : " << P_it << std::endl;
2541 auto P_L2norm_old_axial =
_P_soln->L2norm();
2542 for (
unsigned int iblock = 0; iblock <
_n_blocks; iblock++)
2546 Real T_block_error = 1.0;
2548 _console <<
"Solving Block: " << iblock <<
" From first level: " << first_level
2549 <<
" to last level: " << last_level << std::endl;
2555 _console <<
"Reached maximum number of temperature iterations for block: " << iblock
2559 auto T_L2norm_old_block =
_T_soln->L2norm();
2576 _console <<
"Done with main solve." << std::endl;
2585 _console <<
"Starting thermal solve." << std::endl;
2592 _console <<
"Done with thermal solve." << std::endl;
2597 _console <<
"Start updating thermophysical properties." << std::endl;
2606 _console <<
"Done updating thermophysical properties." << std::endl;
2610 _aux->solution().close();
2612 auto T_L2norm_new =
_T_soln->L2norm();
2614 std::abs((T_L2norm_new - T_L2norm_old_block) / (T_L2norm_old_block + 1E-14));
2615 _console <<
"T_block_error: " << T_block_error << std::endl;
2618 auto P_L2norm_new_axial =
_P_soln->L2norm();
2620 std::abs((P_L2norm_new_axial - P_L2norm_old_axial) / (P_L2norm_old_axial +
_P_out + 1E-14));
2621 _console <<
"P_error :" << P_error << std::endl;
2624 _console <<
"Iteration: " << P_it << std::endl;
2625 _console <<
"Maximum iterations: " << P_it_max << std::endl;
2630 _console <<
"Finished executing subchannel solver\n";
2635 _console <<
"Commencing calculation of Pin surface temperature \n";
2636 for (
unsigned int i_pin = 0; i_pin <
_n_pins; i_pin++)
2638 for (
unsigned int iz = 0; iz <
_n_cells + 1; ++iz)
2642 Real rod_counter = 0.0;
2647 auto mu = (*_mu_soln)(node);
2648 auto S = (*_S_flow_soln)(node);
2649 auto w_perim = (*_w_perim_soln)(node);
2650 auto Dh_i = 4.0 *
S / w_perim;
2651 auto Re = (((*_mdot_soln)(node) /
S) * Dh_i /
mu);
2654 auto Pr = (*_mu_soln)(node)*
cp /
k;
2656 auto hw = Nu *
k / Dh_i;
2658 mooseError(
"Dpin should not be null or negative when computing pin powers: ",
2661 (*_q_prime_soln)(pin_node) / ((*
_Dpin_soln)(pin_node)*M_PI * hw) + (*_T_soln)(node);
2664 if (rod_counter > 0)
2665 _Tpin_soln->set(pin_node, sumTemp / rod_counter);
2667 mooseError(
"Pin was not found for pin index: " + std::to_string(i_pin));
2675 _console <<
"Commencing calculation of duct surface temperature " << std::endl;
2677 for (Node * dn : duct_nodes)
2680 auto mu = (*_mu_soln)(node_chan);
2681 auto S = (*_S_flow_soln)(node_chan);
2682 auto w_perim = (*_w_perim_soln)(node_chan);
2683 auto Dh_i = 4.0 *
S / w_perim;
2684 auto Re = (((*_mdot_soln)(node_chan) /
S) * Dh_i /
mu);
2687 auto Pr = (*_mu_soln)(node_chan)*
cp /
k;
2689 auto hw = Nu *
k / Dh_i;
2695 _aux->solution().close();
2698 Real power_in = 0.0;
2699 Real power_out = 0.0;
2700 Real Total_surface_area = 0.0;
2701 Real Total_wetted_perimeter = 0.0;
2702 Real mass_flow_in = 0.0;
2703 Real mass_flow_out = 0.0;
2704 for (
unsigned int i_ch = 0; i_ch <
_n_channels; i_ch++)
2708 Total_surface_area += (*_S_flow_soln)(node_in);
2709 Total_wetted_perimeter += (*_w_perim_soln)(node_in);
2710 power_in += (*_mdot_soln)(node_in) * (*
_h_soln)(node_in);
2711 power_out += (*_mdot_soln)(node_out) * (*
_h_soln)(node_out);
2712 mass_flow_in += (*_mdot_soln)(node_in);
2713 mass_flow_out += (*_mdot_soln)(node_out);
2715 auto h_bulk_out = power_out / mass_flow_out;
2716 auto T_bulk_out =
_fp->T_from_p_h(
_P_out, h_bulk_out);
2718 Real bulk_Dh = 4.0 * Total_surface_area / Total_wetted_perimeter;
2720 Real bulk_Re = mass_flow_in * bulk_Dh / (inlet_mu * Total_surface_area);
2723 _console <<
" ======================================= " << std::endl;
2724 _console <<
" ======== Subchannel Print Outs ======== " << std::endl;
2725 _console <<
" ======================================= " << std::endl;
2726 _console <<
"Total flow area :" << Total_surface_area <<
" m^2" << std::endl;
2727 _console <<
"Assembly hydraulic diameter :" << bulk_Dh <<
" m" << std::endl;
2728 _console <<
"Assembly Re number :" << bulk_Re <<
" [-]" << std::endl;
2729 _console <<
"Bulk coolant temperature at outlet :" << T_bulk_out <<
" K" << std::endl;
2730 _console <<
"Power added to coolant is : " << power_out - power_in <<
" Watt" << std::endl;
2731 _console <<
"Mass flow rate in is : " << mass_flow_in <<
" kg/sec" << std::endl;
2732 _console <<
"Mass balance is : " << mass_flow_out - mass_flow_in <<
" kg/sec" << std::endl;
2733 _console <<
"User defined outlet pressure is : " <<
_P_out <<
" Pa" << std::endl;
2734 _console <<
" ======================================= " << std::endl;
2739 "Energy conservation equation might not be solved correctly, Power added to coolant: " +
2740 std::to_string(power_out - power_in) +
" Watt ");
Vec _amc_friction_force_rhs
Mat _amc_sys_mdot_mat
Axial momentum system matrix.
static const std::string PRESSURE_DROP
pressure drop
static const std::string MASS_FLOW_RATE
mass flow rate
const bool _pin_mesh_exist
Flag that informs if there is a pin mesh or not.
void computeRho(int iblock)
Computes Density per channel for block iblock.
Vec _amc_gravity_rhs
Axial momentum conservation - buoyancy force No implicit matrix.
virtual const unsigned int & getNumOfChannels() const =0
Return the number of channels per layer.
std::unique_ptr< SolutionHandle > _Tduct_soln
static const std::string DENSITY
density
void computeSumWij(int iblock)
Computes net diversion crossflow per channel for block iblock.
unsigned int _n_blocks
number of axial blocks
const bool _compute_power
Flag that informs if we need to solve the Enthalpy/Temperature equations or not.
virtual const Real & getPinDiameter() const
Return Pin diameter.
static InputParameters validParams()
std::unique_ptr< SolutionHandle > _T_soln
std::unique_ptr< SolutionHandle > _h_soln
const PetscReal & _dtol
The divergence tolerance for the ksp linear solver.
virtual Real computeBeta(unsigned int i_gap, unsigned int iz)=0
Computes turbulent mixing coefficient.
const bool _monolithic_thermal_bool
Thermal monolithic bool.
std::unique_ptr< SolutionHandle > _q_prime_duct_soln
virtual const std::vector< unsigned int > & getPinChannels(unsigned int i_pin) const =0
Return a vector of channel indices for a given Pin index.
virtual Node * getChannelNodeFromDuct(Node *channel_node)=0
Function that gets the channel node from the duct node.
T & getMesh(MooseMesh &mesh)
function to cast mesh
PetscScalar _max_sumWij_new
Mat _cmc_friction_force_mat
Cross momentum conservation - friction force.
void computeP(int iblock)
Computes Pressure per channel for block iblock.
Real _TR
Flag that activates or deactivates the transient parts of the equations we solve by multiplication...
libMesh::DenseMatrix< Real > _Wij_old
const PostprocessorValue & _P_out
Outlet Pressure.
SubChannelMesh & _subchannel_mesh
void computeT(int iblock)
Computes Temperature per channel for block iblock.
virtual void zero() override final
virtual const std::vector< Real > & getZGrid() const
Get axial location of layers.
Vec _amc_pressure_force_rhs
void computeMdot(int iblock)
Computes mass flow per channel for block iblock.
virtual const Real & getPitch() const
Return the pitch between 2 subchannels.
virtual EChannelType getSubchannelType(unsigned int index) const =0
Return the type of the subchannel for given subchannel index.
std::unique_ptr< SolutionHandle > _S_flow_soln
Vec _hc_time_derivative_rhs
const bool _compute_density
Flag that activates or deactivates the calculation of density.
Vec _mc_axial_convection_rhs
Mat _amc_friction_force_mat
Axial momentum conservation - friction force.
static const std::string PIN_DIAMETER
pin diameter
static const std::string DUCT_LINEAR_HEAT_RATE
duct linear heat rate
PetscScalar _added_K
Added resistances for monolithic convergence.
PetscErrorCode createPetscVector(Vec &v, PetscInt n)
Petsc Functions.
static InputParameters validParams()
Vec _cmc_advective_derivative_rhs
const bool _duct_mesh_exist
Flag that informs if there is a pin mesh or not.
static InputParameters validParams()
const Real & _T_tol
Convergence tolerance for the temperature loop in external solve.
Vec _amc_time_derivative_rhs
void computeMu(int iblock)
Computes Viscosity per channel for block iblock.
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
virtual Node * getPinNode(unsigned int i_pin, unsigned iz) const =0
Get the pin mesh node for a given pin index and elevation index.
PetscScalar computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet=0.0)
bool isRestarting() const
const bool _segregated_bool
Segregated solve.
Mat _cmc_sys_Wij_mat
Lateral momentum system matrix.
Vec _amc_advective_derivative_rhs
bool _converged
Variable that informs whether we exited external solve with a converged solution or not...
virtual const std::string & name() const
virtual const std::vector< unsigned int > & getChannelGaps(unsigned int i_chan) const =0
Return a vector of gap indices for a given channel index.
Vec _hc_cross_derivative_rhs
const bool _staggered_pressure_bool
Flag to define the usage of staggered or collocated pressure.
virtual void initializeSolution()=0
Function to initialize the solution & geometry fields.
void mooseWarning(Args &&... args) const
virtual const MooseVariableFieldBase & getVariable(const THREAD_ID tid, const std::string &var_name, Moose::VarKindType expected_var_type=Moose::VarKindType::VAR_ANY, Moose::VarFieldType expected_var_field_type=Moose::VarFieldType::VAR_FIELD_ANY) const override
Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz)
Function that computes the heat flux added by the duct.
Mat _mc_sumWij_mat
Matrices and vectors to be used in implicit assembly Mass conservation Mass conservation - sum of cro...
PetscErrorCode populateVectorFromDense(Vec &x, const T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
static const std::string DUCT_TEMPERATURE
duct temperature
SubChannel1PhaseProblem * schp
const Real & _CT
Turbulent modeling parameter used in axial momentum equation.
virtual Node * getChannelNode(unsigned int i_chan, unsigned iz) const =0
Get the subchannel mesh node for a given channel index and elevation index.
static const std::string WETTED_PERIMETER
wetted perimeter
virtual void computeh(int iblock)=0
Computes Enthalpy per channel for block iblock.
std::unique_ptr< SolutionHandle > _rho_soln
Mat _amc_turbulent_cross_flows_mat
Mass conservation - density time derivative No implicit matrix.
static const std::string cp
Mat _amc_advective_derivative_mat
Axial momentum conservation - advective (Eulerian) derivative.
static const std::string VISCOSITY
viscosity
Vec _hc_added_heat_rhs
Enthalpy conservation - source and sink.
processor_id_type n_processors() const
std::vector< Real > _z_grid
axial location of nodes
static const std::string PIN_TEMPERATURE
pin temperature
Vec _cmc_time_derivative_rhs
const int & _T_maxit
Maximum iterations for the inner temperature loop.
Mat _amc_pressure_force_mat
Axial momentum conservation - pressure force.
static const std::string LINEAR_HEAT_RATE
linear heat rate
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
void computeWijResidual(int iblock)
Computes Residual Matrix based on the lateral momentum conservation equation for block iblock...
const std::vector< double > x
static const std::string S
Real f(Real x)
Test function for Brents method.
std::unique_ptr< SolutionHandle > _q_prime_soln
virtual void syncSolutions(Direction direction) override
static const std::string mu
static const std::string pitch
const PetscReal & _atol
The absolute convergence tolerance for the ksp linear solver.
PetscScalar computeInterpolationCoefficients(PetscScalar Peclet=0.0)
Functions that computes the interpolation scheme given the Peclet number.
static const std::string ENTHALPY
enthalpy
Real root(std::function< Real(Real)> const &f, Real x1, Real x2, Real tol=1.0e-12)
Finds the root of a function using Brent's method.
std::shared_ptr< AuxiliarySystem > _aux
virtual const unsigned int & getNumOfPins() const =0
Return the number of pins.
Mat _hc_cross_derivative_mat
Enthalpy conservation - cross flux derivative.
Mat _mc_axial_convection_mat
Mass conservation - axial convection.
virtual Node * getDuctNodeFromChannel(Node *channel_node)=0
Function that gets the duct node from the channel node.
void initialSetup() override
auto Peclet(const T1 &volume_fraction, const T2 &cp, const T3 &rho, const T4 &vel, const T5 &D_h, const T6 &k)
Compute Peclet number.
virtual Real computeFrictionFactor(FrictionStruct friction_args)=0
Returns friction factor.
virtual ~SubChannel1PhaseProblem()
const PetscInt & _maxit
The maximum number of iterations to use for the ksp linear solver.
PetscErrorCode createPetscMatrix(Mat &M, PetscInt n, PetscInt m)
friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void *ctx)
This is the residual Vector function in a form compatible with the SNES PETC solvers.
virtual const std::vector< std::vector< Real > > & getKGrid() const
Get axial cell location and value of loss coefficient.
const MooseEnum _interpolation_scheme
The interpolation method used in constructing the systems.
std::unique_ptr< SolutionHandle > _Dpin_soln
Base class for the 1-phase steady-state/transient subchannel solver.
libMesh::DenseMatrix< Real > _Wij_residual_matrix
std::unique_ptr< SolutionHandle > _mdot_soln
PetscErrorCode populateDenseFromVector(const Vec &x, T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
static const std::string SUM_CROSSFLOW
sum of diversion crossflow
virtual const unsigned int & getNumOfAxialCells() const
Return the number of axial cells.
virtual Real getGapWidth(unsigned int axial_index, unsigned int gap_index) const =0
Return gap width for a given gap index.
LibmeshPetscCallQ(DMShellGetContext(dm, &ctx))
const PetscReal & _rtol
The relative convergence tolerance, (relative decrease) for the ksp linear solver.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unique_ptr< SolutionHandle > _displacement_soln
PetscErrorCode formFunction(SNES, Vec x, Vec f, void *ctx)
virtual const std::vector< Node * > getDuctNodes() const =0
Function that return the vector with the maps to the nodes if they exist.
std::unique_ptr< SolutionHandle > _DP_soln
const bool _compute_viscosity
Flag that activates or deactivates the calculation of viscosity.
static const std::string PRESSURE
pressure
libMesh::DenseMatrix< Real > & _Wij
virtual void externalSolve() override
void computeWijPrime(int iblock)
Computes turbulent crossflow per gap for block iblock.
static const std::string alpha
virtual const Real & getCrossflowSign(unsigned int i_chan, unsigned int i_local) const =0
Return a sign for the crossflow given a subchannel index and local neighbor index.
std::unique_ptr< SolutionHandle > _SumWij_soln
PetscErrorCode petscSnesSolver(int iblock, const libMesh::DenseVector< Real > &solution, libMesh::DenseVector< Real > &root)
Computes solution of nonlinear equation using snes and provided a residual in a formFunction.
Vec _hc_advective_derivative_rhs
libMesh::DenseVector< Real > residualFunction(int iblock, libMesh::DenseVector< Real > solution)
Computes Residual Vector based on the lateral momentum conservation equation for block iblock & updat...
static const std::string SURFACE_AREA
surface area
bool absoluteFuzzyLessEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
void resize(const unsigned int new_m, const unsigned int new_n)
const bool _verbose_subchannel
Boolean to printout information related to subchannel solve.
PetscScalar _correction_factor
Mat _amc_cross_derivative_mat
Axial momentum conservation - cross flux derivative.
const Real & _P_tol
Convergence tolerance for the pressure loop in external solve.
void mooseError(Args &&... args) const
virtual bool solverSystemConverged(const unsigned int) override
virtual const unsigned int & getNumOfGapsPerLayer() const =0
Return the number of gaps per layer.
struct SubChannel1PhaseProblem::FrictionStruct _friction_args
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
Mat _hc_advective_derivative_mat
Enthalpy conservation - advective (Eulerian) derivative;.
Base class for subchannel meshes.
void computeWijFromSolve(int iblock)
Computes diversion crossflow per gap for block iblock.
const bool _implicit_bool
Flag to define the usage of a implicit or explicit solution.
Mat _amc_time_derivative_mat
Axial momentum conservation - time derivative.
const ConsoleStream _console
static const std::string DISPLACEMENT
subchannel displacement
Mat _hc_time_derivative_mat
Enthalpy Enthalpy conservation - time derivative.
virtual bool isTransient() const override
PetscErrorCode implicitPetscSolve(int iblock)
Computes implicit solve using PetSc.
void computeDP(int iblock)
Computes Pressure Drop per channel for block iblock.
SubChannel1PhaseProblem(const InputParameters ¶ms)
Mat _cmc_time_derivative_mat
Cross momentum Cross momentum conservation - time derivative.
libMesh::DenseMatrix< Real > _WijPrime
virtual void initialSetup() override
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
Vec _cmc_friction_force_rhs
bool isRecovering() const
static const std::string TEMPERATURE
temperature
Mat _cmc_pressure_force_mat
Cross momentum conservation - pressure force.
virtual Real & dt() const
Vec _amc_turbulent_cross_flows_rhs
MooseUnits pow(const MooseUnits &, int)
std::unique_ptr< SolutionHandle > _P_soln
static const std::string k
std::unique_ptr< SolutionHandle > _w_perim_soln
void ErrorVector unsigned int
Vec _cmc_pressure_force_rhs
Vec _amc_cross_derivative_rhs
Mat _cmc_advective_derivative_mat
Cross momentum conservation - advective (Eulerian) derivative.
const SinglePhaseFluidProperties * _fp
Solutions handles and link to TH tables properties.
Mat _hc_sys_h_mat
System matrices.
std::unique_ptr< SolutionHandle > _mu_soln
virtual const std::pair< unsigned int, unsigned int > & getGapChannels(unsigned int i_gap) const =0
Return a pair of subchannel indices for a given gap index.
std::unique_ptr< SolutionHandle > _Tpin_soln