Some incomprehensible parts of SCM code

[Rforoeverad]

$\qquad$ When I read the SCM source code, I encounter the following difficult problems to understand. I hope that developers can take time to give a reply.

1. When the implicit segregated method is used, in the SubChannel1PhaseProblem::computeWijResidual(int iblock) member function, the member variable _Wij_vec is always a 0 vector and loc_holder_Wij is destroyed before it is used.
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 1794 to 1796 in 832d6f3

   	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
   	loc_holder_Wij, _Wij, first_node, last_node, _n_gaps));
   	LibmeshPetscCall(MatMult(_cmc_sys_Wij_mat, _Wij_vec, sol_holder_W));

2. When the implicit monolithic method is used, in the SubChannel1PhaseProblem::implicitPetscSolve(int iblock) member function, the bool variable lag_block_thermal_solve is always ture such that loc_vecs[3] is always solved as a 1 vector. The result is logically the same as the case where _monolithic_thermal_bool is false. So why not directly solve _h_soln in the global matrix?
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 2095 to 2111 in 832d6f3

   	// Energy conservation
   	if (_monolithic_thermal_bool)
   	{
   	field_num = 3;
   	mat_array[Q * field_num + 0] = NULL;
   	mat_array[Q * field_num + 1] = NULL;
   	mat_array[Q * field_num + 2] = NULL;
   	LibmeshPetscCall(MatDuplicate(_hc_sys_h_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 3]));
   	if (lag_block_thermal_solve)
   	{
   	LibmeshPetscCall(MatZeroEntries(mat_array[Q * field_num + 3]));
   	LibmeshPetscCall(MatShift(mat_array[Q * field_num + 3], 1.0));
   	}
   	LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
   	LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
   	LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &vec_array[field_num]));
   	LibmeshPetscCall(VecCopy(_hc_sys_h_rhs, vec_array[field_num]));

3. Why does SCM use one-dimensional mesh to solve, which is consistent with the principle of CTF code?
   R.Salko et al. CTF Theory Manual: Version 4.4. 2023

The convergence eficiency can also be greatly enhanced with the use of rebalancing. This is the process of reducing the mesh to a 1D mesh and only solving pressure changes in the axial direction. Direct inversion is used. but the number of mesh cells is greatly reduced because the only concern is the axial direction.

$\qquad$ Many thanks for your help!!

[GiudGiud]

@kyriv1980

[kyriv1980]

Regarding the first question : _Wij_vec is the current solution of the crossflow. Its a PETSc vector that has the size of all axial levels times gaps per level. It's not a zero vector. loc_holder_Wij does seem not to be used. So that might be a bug. @tanoret

[Rforoeverad]

$\qquad$ _Wij_vec is updated only when it is used to relax Wij in implicitPetscSolve, then _Wij_vec is never updated when using the implicit segregated method, so it maintains the zero vector value created using the createPetscVector function. If so, even if _Wij_vec represents the current solution of the crossflow, it can be temporarily defined in implicitPetscSolve without being a class member variable.

[Rforoeverad]

$\qquad$ In the case of variables defined but not used, is the variable _cmc_Wij_channel_dummy defined in the SubChannel1PhaseProblem.h file created to define the boundary conditions and what is the role of the variable _global_counter and _dpz_error?

[kyriv1980]

computeWijFromSolve(iblock) computes _Wij for the segregated methods (implicit-segregated, explicit-segregated). _Wij_vec is an intermediate way to store _Wij in a type that is usable in the implicit solve. It probably doesn't need to be a member variable (_Wij_vec) but it works.

[kyriv1980]

The purpose of _cmc_Wij_channel_dummy has been lost in time. I should probably delete it. Thanks for catching that. The same applies for _global_counter. I can't find _dpz_error ... maybe you'll have to rebase on devel???

[Rforoeverad]

The location of _dpz_error is shown here. Similarly, the structure StructPetscMatVec is not used.
InterWrapper1PhaseProblem may have the same problems mentioned above.

moose/modules/subchannel/include/problems/SubChannel1PhaseProblem.h

Lines 115 to 116 in 832d6f3

	/// average relative error in pressure drop of channels
	Real _dpz_error;

moose/modules/subchannel/include/problems/SubChannel1PhaseProblem.h

Lines 40 to 44 in 832d6f3

	struct StructPetscMatVec
	{
	Mat A;
	Vec x;
	};

[kyriv1980]

Regarding the second question. It seems you are right regarding the enthalpy calculation. It is always lagged. @tanoret ?

[Rforoeverad]

$\qquad$ The similar question is that the three bool variables _axial_mass_flow_tight_coupling, _pressure_axial_momentum_tight_coupling and _pressure_cross_momentum_tight_coupling used to characterize coupling are always defined as true in the SubChannel1PhaseProblem.C file. So their role is only to indicate the presence of coupling effect, or say there are other development purposes?

[kyriv1980]

I think that the purpose of the developer here was to make the code expandable and versatile in future iterations.

[Rforoeverad]

OK, thanks for your response.

[kyriv1980]

Regarding the third question: SCM uses 1D mesh per sub-channel that is radially connected to the neighboring 1D channels, through the lateral momentum equation to capture the 3D effect of crossflow.

[Rforoeverad]

$\qquad$ Does the selection of Wij as a JFNK nonlinear solution variable also take into account the introduction of strong coupling characteristics due to the need to describe the three-dimensional characteristics of the variable?

[kyriv1980]

The selection of _Wij, stems intuitively. Given _Wij the axial mass flow can be retrieved for all nodes. With the axial mass flow known pressure drop and pressure can be retrieved.

Then the process can repeat.

[Rforoeverad]

Thanks a lot.

[Rforoeverad]

$\qquad$ There are some differences between the formula and the code when calculating the specific form of the following variables:

1. I think that the former area of mass_term1 of computeDP should be (*_S_flow_soln)(node_out).
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 692 to 705 in 832d6f3

   	auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
   	auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
   	auto rho_in = (*_rho_soln)(node_in);
   	auto rho_out = (*_rho_soln)(node_out);
   	auto mu_in = (*_mu_soln)(node_in);
   	auto S = (*_S_flow_soln)(node_in);
   	auto w_perim = (*_w_perim_soln)(node_in);
   	// hydraulic diameter in the i direction
   	auto Dh_i = 4.0 * S / w_perim;
   	auto time_term = _TR * ((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
   	dz * 2.0 * (*_mdot_soln)(node_out) * (rho_out - _rho_soln->old(node_out)) /
   	rho_in / _dt;
   	auto mass_term1 =
   	std::pow((*_mdot_soln)(node_out), 2.0) * (1.0 / S / rho_out - 1.0 / S / rho_in);

2. The mass_term_out and mass_term_in of the inertia_term of the function computeWijResidual are obtained from the same axial level. Does this mean that the use of the axial area of the inlet and outlet has little effect on the results?
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 1438 to 1444 in 8b12082

   	auto mass_term_out =
   	(*_mdot_soln)(node_out_i) / Si / rho_i + (*_mdot_soln)(node_out_j) / Sj / rho_j;
   	auto mass_term_in =
   	(*_mdot_soln)(node_in_i) / Si / rho_i + (*_mdot_soln)(node_in_j) / Sj / rho_j;
   	auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out * _Wij(i_gap, iz);
   	auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in * _Wij(i_gap, iz - 1);
   	auto inertia_term = term_out - term_in;

3. The thermal conductivity of e_cond of the function computeh uses the arithmetic mean. Will it be more in line with the physical reality to consider using the harmonic mean?
   

   moose/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C

   Lines 788 to 800 in 832d6f3

   	auto shape_factor =
   	0.66 * (pitch / pin_diameter) *
   	std::pow((_subchannel_mesh.getGapWidth(iz, i_gap) / pin_diameter), -0.3);
   	if (ii_ch == i_ch)
   	{
   	e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
   	((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
   	}
   	else
   	{
   	e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
   	((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
   	}

4. Why does sweep_enthalpy in computeh need to be listed separately? It can be completely listed in the Cheng-Todreas model.
   

   moose/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C

   Lines 734 to 741 in 832d6f3

   	// Calculation of turbulent mixing parameter
   	auto beta = Cs * std::pow(Ar2 / A2, 0.5) * std::tan(theta);
   	auto wsweep_in = gedge_ave * beta * Sij;
   	auto wsweep_out = gedge_ave * beta * Sij;
   	auto sweep_hin = (*_h_soln)(node_sin);
   	auto sweep_hout = (*_h_soln)(node_in);
   	sweep_enthalpy = (wsweep_in * sweep_hin - wsweep_out * sweep_hout);

5. Compared with the calculation of other upwind values, the h_star variable of the computeh function does not consider the case of _Wij = 0, so why write this case in particular?
   

   moose/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C

   Line 1153 in 832d6f3

   else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary

   
   $\qquad$ Thank you very much for your reply during your busy schedule!!

[GiudGiud]

The thermal conductivity of e_cond of the function computeh uses the arithmetic mean. Will it be more in line with the physical reality to consider using the harmonic mean ?

Yes but here the differences are expected to be small. Subchannel simulations are ran with a fairly fine axial mesh. It is not like the porous media solver in Navier Stokes where different porosity-regions might have different thermal conductivities.

[kyriv1980]

• When first developing SCM, subchannel surface area was thought to be constant so the node selections doesn't really affect the results for that case. In any case we usually opt to upwind surface area which is consistent with the one-directional flow assumption in most subchannel codes.

• Maybe in future iterations we can have a more consistent indexing scheme for the various geometric parameters. If you want, correct it where you see fit and make a PR about it.

• The sweep mixing parameter is only considered for enthalpy and not momentum, which is not the case for the SCM where in general the mixing parameter beta affects both momentum mixing and enthalpy mixing. So sweep flow had to be separated and the momentum mixing set to zero.... but we plan to refactor in the future so that the implementation becomes more consistent with the rest of the code.

• if _Wij = zero, all the terms that we would add to the matrices/vectors become zero so there is no need to consider this option.

[Rforoeverad]

The part about mass_term1 appears in this dissertation.

[kyriv1980]

The equation shown here is used for the explicit method only. The implicit solves discretize the inertia term differently.

[Rforoeverad]

Many thanks!!

[Rforoeverad]

$\qquad$ There are some questions about the interpolation as follows.
$\qquad$ When calculating the transverse momentum equation, the control volume used is staggered by half the grid with the axial control volume in the case of _staggered_pressure_bool assigned to false, and the interpolation is related to the spatial discretization and has nothing to do with the explicit and implicit processing. So why does the implicit case need interpolation and the explicit case not use interpolation? The same question appears in computeDP, computeP and computeWijResidual.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 1414 to 1425 in 8b12082

	auto chans = _subchannel_mesh.getGapChannels(i_gap);
	unsigned int i_ch = chans.first;
	unsigned int j_ch = chans.second;
	auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
	auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
	auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
	auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
	auto rho_i = (*_rho_soln)(node_in_i);
	auto rho_j = (*_rho_soln)(node_in_j);
	auto Si = (*_S_flow_soln)(node_in_i);
	auto Sj = (*_S_flow_soln)(node_in_j);
	auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 1484 to 1496 in 8b12082

	// inlet, outlet, and interpolated densities
	auto rho_i_in = (*_rho_soln)(node_in_i);
	auto rho_i_out = (*_rho_soln)(node_out_i);
	auto rho_i_interp = computeInterpolatedValue(rho_i_out, rho_i_in, 0.5);
	auto rho_j_in = (*_rho_soln)(node_in_j);
	auto rho_j_out = (*_rho_soln)(node_out_j);
	auto rho_j_interp = computeInterpolatedValue(rho_j_out, rho_j_in, 0.5);
	// inlet, outlet, and interpolated areas
	auto S_i_in = (*_S_flow_soln)(node_in_i);
	auto S_i_out = (*_S_flow_soln)(node_out_i);
	auto S_j_in = (*_S_flow_soln)(node_in_j);
	auto S_j_out = (*_S_flow_soln)(node_out_j);

[kyriv1980]

The explicit was the first method developed. Since it assumes one directional flow, most needed quantities are up-winded.

[Rforoeverad]

I'm deeply grateful for your assistance!!

[Rforoeverad]

$\qquad$ There may be a small problem in the nonlinear solution part.
$\qquad$ In the petscSnesSolver function, when Vec x is created but not initialized, x is initialized to 0 vector by default in my machine, and the default initialization value may generate different default initialization values on different machines. Then, the formFunction(SNES, Vec x, Vec f, void * ctx) requires x as a function parameter. Will this cause problems caused by default initialization parameters?

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 1883 to 1897 in 8b12082

	Ctx ctx;
	ctx.iblock = iblock;
	ctx.schp = this;
	LibmeshPetscCall(SNESSetFunction(snes, r, formFunction, &ctx));
	LibmeshPetscCall(SNESGetKSP(snes, &ksp));
	LibmeshPetscCall(KSPGetPC(ksp, &pc));
	LibmeshPetscCall(PCSetType(pc, PCNONE));
	LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
	LibmeshPetscCall(SNESSetFromOptions(snes));
	LibmeshPetscCall(VecGetArray(x, &xx));
	for (unsigned int i = 0; i < _block_size * _n_gaps; i++)
	{
	xx[i] = solution(i);
	}
	LibmeshPetscCall(VecRestoreArray(x, &xx));

[kyriv1980]

In the PETSc user manual that i used as a reference to built the SNES solver in SCM, the x vector which is the initial guess/current iterate, gets populated just before the solve, which is what i do here : https://petsc.org/release/manual/snes/#sec-nlmatrixfree

[Rforoeverad]

OK, this may be my cognitive deficit. Thank you again for taking the time to give a patient and professional reply!!

[Rforoeverad]

$\qquad$ There are some questions about Wij as follows.

1. According to the annotations and related calculation logic, Wij is defined as positive from the small numbered subchannel to the large numbered subchannel. So why not define the direction according to the actual physical situation, that is, the orientation of the coordinate system, so that the transverse flow in the actual physical situation can be obtained?
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 483 in 8b12082

   // Wij positive out of i into j;

   
   

   moose/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C

   Line 750 in 832d6f3

   // i is always the smallest and first index in the mapping

   
   

   moose/modules/subchannel/include/mesh/TriSubChannelMesh.h

   Lines 195 to 199 in 8b12082

   	/// Defines the global cross-flow direction -1 or 1 for each subchannel and
   	/// for all gaps that are belonging to the corresponding subchannel.
   	/// Given a subchannel and a gap, if the neighbor subchannel index belonging to the same gap is lower,
   	/// set it to -1, otherwise set it to 1.
   	std::vector<std::vector<Real>> _sign_id_crossflow_map;

2. In the computeDP and computeh functions, cross_index is used to determine the upwind value in some places, while the value of Wij at iz is used elsewhere. Maybe it is possible to unify the description, or does it show that it is more accurate to use the existing or interpolated Wij to determine the upwind direction at the same axial position as *_star?
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 623 in 8b12082

   if (_Wij(i_gap, iz) > 0.0)

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 798 in 8b12082

   unsigned int cross_index = iz; // iz-1;

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 818 in 8b12082

   if (_Wij(i_gap, cross_index) > 0.0)

3. sol_Wij_array and sumWij_loc in implicitPetscSolve may not be required. In addition, the function VecRestoreArray may be required after the variable sol_p_array and sol_h_array are used.
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 2388 to 2389 in 8b12082

   	PetscScalar * sol_Wij_array;
   	LibmeshPetscCall(VecGetArray(sol_Wij, &sol_Wij_array));

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 2484 to 2487 in 8b12082

   	Vec sumWij_loc;
   	LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels));
   	LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
   	_prod, *_SumWij_soln, first_node, last_node, _n_channels));

[kyriv1980]

Actually the comment: "i is the smallest and first index in the matrix" is false. _Wij is considered positive if it's coming out of channel i and negative otherwise.

Cross_index and iz are usually interchangeable.

[Rforoeverad]

Okay, I'll dig a little bit deeper into the code later, and finally thanks for providing such an elegant code!!

[Rforoeverad]

Hello, please forgive me for bothering you again.
I write the following changes here and hope you can take the time to look over:

1. Here are some changes for the populate function template.
   When using the populate template function, some places are inconsistent when updating variables.
   

   moose/modules/subchannel/include/problems/SubChannel1PhaseProblem.h

   Line 382 in 8b12082

   for (unsigned int iz = first_axial_level; iz < last_axial_level; iz++)

for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)

$\qquad$ This modification is to keep the value range of the populateVectorFromDense template function consistent with other populate template function.
$\qquad$ On the basis of the above modification, I modify the following populate function call statement. Please check whether the relevant calls of populate function achieve the desired purpose.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 456 to 457 in 8b12082

	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
	loc_Wij, _Wij, first_node, last_node + 1, _n_gaps));

LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
          loc_Wij, _Wij, first_node, last_node, _n_gaps));

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 2408 to 2409 in 8b12082

	LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
	sol_Wij, _Wij, first_node, last_node - 1, _n_gaps));

LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
      sol_Wij, _Wij, first_node, last_node, _n_gaps));

$\qquad$ In addition, populateSolutionGap is not used and the remaining four populate template functions can be written as member functions. After all, the template parameter type T can be inferred from calls and names.

2. This is a related modification to the previously mentioned _Wij_vec.
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 1697 in 8b12082

   LibmeshPetscCall(MatMult(_cmc_sys_Wij_mat, _Wij_vec, sol_holder_W));

LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
          _Wij_vec, _Wij, first_node + 1, last_node, _n_gaps));
LibmeshPetscCall(MatMult(_cmc_sys_Wij_mat, _Wij_vec, sol_holder_W));

$\qquad$ As previously discussed, in the implicit segregated calculations, _Wij_vec needs to be updated here logically, otherwise it is always a 0 vector value.
$\qquad$ Note that the first_node in this member function is iblock * _block_size and the range of value remains associated with the populateVectorFromDense template function defined.

3. Here is the modification of some explicit calculations.
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 602 to 603 in 8b12082

   	auto mass_term1 =
   	std::pow((*_mdot_soln)(node_out), 2.0) * (1.0 / S / rho_out - 1.0 / S / rho_in);

auto mass_term1 = std::pow((*_mdot_soln)(node_out), 2.0) *
                          (1.0 / (*_S_flow_soln)(node_out) / rho_out - 1.0 / S / rho_in);

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 1438 to 1439 in 8b12082

	auto mass_term_out =
	(*_mdot_soln)(node_out_i) / Si / rho_i + (*_mdot_soln)(node_out_j) / Sj / rho_j;

auto mass_term_out =
            (*_mdot_soln)(node_out_i) / (*_S_flow_soln)(node_out_i) / (*_rho_soln)(node_out_i) +
            (*_mdot_soln)(node_out_j) / (*_S_flow_soln)(node_out_j) / (*_rho_soln)(node_out_j);

$\qquad$ The mass_term1 and inertia_term may need to be modified according to the discretization of the axial gradient term in the dissertation mentioned above. But There are two versions of axial velocity $u_{i,k}$ calculation in the code. This may be the index position that needs to be standardized.
$u_{i,k} = \frac{m_{i,k}}{\rho_{i,k}*S_{i,k-1}}$
$u_{i,k} = \frac{m_{i,k}}{\rho_{i,k-1}*S_{i,k-1}}$

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 599 to 601 in 8b12082

	auto time_term = _TR * ((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
	dz * 2.0 * (*_mdot_soln)(node_out) * (rho_out - _rho_soln->old(node_out)) /
	rho_in / _dt;

auto time_term = _TR * (((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
                                dz * 2.0 * (*_mdot_soln)(node_out) *
                                    (rho_out - _rho_soln->old(node_out)) / rho_in / _dt);

$\qquad$ The second term in time_term may also be zero at steady state, so the parentheses after _TR contain the entire line of code.

4. Here are some annotations that may need to be modified.
   

   moose/modules/subchannel/include/problems/SubChannel1PhaseProblem.h

   Lines 125 to 126 in 8b12082

   	/// Flag that informs if there is a pin mesh or not
   	const bool _duct_mesh_exist;

   
   

   moose/modules/subchannel/include/problems/SubChannel1PhaseProblem.h

   Lines 137 to 138 in 8b12082

   	/// Convergence tolerance for the temperature loop in external solve
   	const Real & _T_tol;

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Lines 76 to 77 in 8b12082

   	params.addParam<bool>(
   	"staggered_pressure", false, "Boolean to define the use of explicit or implicit solution.");

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 428 in 8b12082

   // Add to matrix if explicit

   
   

   moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

   Line 2643 in 8b12082

   // Calculate sum of pin surface temperatures that the channels around the pin see

Thanks and best wishes.

[Rforoeverad]

$\qquad$ Thank you for your detailed reply, which has benefited me!! But there are still some inconsistencies.@kyriv1980
1.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 391 to 399 in 5dfaae1

	int i = 0;
	for (unsigned int iz = first_node; iz < last_node + 1; iz++)
	{
	for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
	{
	_Wij(i_gap, iz) = root(i);
	i++;
	}
	}

$\qquad$ The segregated method updates _Wij from first_node to last_node.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 2408 to 2409 in 5dfaae1

	LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
	sol_Wij, _Wij, first_node, last_node - 1, _n_gaps));

$\qquad$ The monolithic method updates _Wij from first_node to last_node - 1.
2.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 456 to 457 in 5dfaae1

	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
	loc_Wij, _Wij, first_node, last_node + 1, _n_gaps));

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 2192 to 2196 in 5dfaae1

	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
	_Wij_loc_vec, _Wij, first_node, last_node, _n_gaps));
	LibmeshPetscCall(VecAbs(_Wij_loc_vec));
	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
	_Wij_old_loc_vec, _Wij_old, first_node, last_node, _n_gaps));

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 2302 to 2303 in 5dfaae1

	LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
	_Wij_vec, _Wij, first_node, last_node, _n_gaps));

3.

moose/modules/subchannel/src/problems/SubChannel1PhaseProblem.C

Lines 2286 to 2287 in 5dfaae1

	LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
	_prod, *_P_soln, first_node, last_node, _n_channels));

LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
        _prod, *_P_soln, first_node - 1, last_node - 1, _n_channels));

$\qquad$ According to the reversibility of the update process, it may be necessary to consider updating _prod from first_node - 1 to last_node - 1.
$\qquad$ The rest of the modifications may not be worth mentioning an issue. Perhaps when I'm more familiar with SCM or find valuable issues, I will summarize and submit these mentioned above.

[kyriv1980]

-1
there is this loop in the method: populateDenseFromVector

for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)

that's why there is the LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix>(
sol_Wij, _Wij, first_node, last_node - 1, _n_gaps)); in the monolithic solver.

-2
populateVectorFromDense has

for (unsigned int iz = first_axial_level; iz < last_axial_level; iz++)

so all calls should have: last node + 1

In general the matrix _Wij.resize(_n_gaps, _n_cells + 1); Is bigger than it needs to. for axial level 0 it holds values that never get populated or are involved in the solve process. On the other hand the vector version of crossflow has exact sizes. So there is some play with the axial indexes to figure out the correct size.

There might be inconsistencies. If you have time please find them and fix them in a PR.

[kyriv1980]

#31275

[kyriv1980]

I created an issue you can reference if you want to create a PR. Thank you.

[Rforoeverad]

$\qquad$ Thank you for your reply.@kyriv1980
$\qquad$ However, for the first problem, n_blocks is usually 1 in most cases but not all, so when using the monolithic method, the original statement never updates the value at the last_node, that is, outlet, which may lead to non-conservation and is not consistent with the segregated method that updates the value at the last_node.

[Rforoeverad]

$\qquad$ In addition, my graduation thesis design is closely related to SCM, so I'd like to know whether you will develop a newton method for implicit monolithic solution, apply SCM with Inter-wrapper flow models(IWF) to the whole reactor simulation and consider SCM as a component of other open source modules or code recently.
$\qquad$ Any advice will be greatly appreciated!!

[kyriv1980]

1. Yes

2. SCM can be coupled to other moose APPS to create multi-assemly models that have SCM run in parallel. An example case will be uploaded to the VTB soon. https://github.com/idaholab/virtual_test_bed

[Rforoeverad]

Thank you for your information. After that, I will follow up SCM related issues, discussion, history, etc.
Best Regards!!

[kyriv1980]

#31276

[Rforoeverad]

$\qquad$ I have submitted all the modifications(#31283) involved in this discussion on the basis of #31276. I hope you can take time to review it.
$\qquad$ In addition, similar modifications can be made to InterWrapper1PhaseProblem since it is similar to the code statement of SubChannel1PhaseProblem.