diff --git a/modules/combined/test/tests/subchannel_thm_coupling/subchannel.i b/modules/combined/test/tests/subchannel_thm_coupling/subchannel.i index 20f46a2d3f60..ee33d6ddad8d 100644 --- a/modules/combined/test/tests/subchannel_thm_coupling/subchannel.i +++ b/modules/combined/test/tests/subchannel_thm_coupling/subchannel.i @@ -105,7 +105,6 @@ heated_length = 1.0 implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_multiapps = true verbose_subchannel = false interpolation_scheme = 'upwind' diff --git a/modules/subchannel/doc/content/modules/subchannel/general/subchannel_theory.md b/modules/subchannel/doc/content/modules/subchannel/general/subchannel_theory.md index f052e5235e34..f44ee32ee3e1 100644 --- a/modules/subchannel/doc/content/modules/subchannel/general/subchannel_theory.md +++ b/modules/subchannel/doc/content/modules/subchannel/general/subchannel_theory.md @@ -118,7 +118,7 @@ h_{ij}' = \sum_{j} w_{ij}'\Delta h_{ij} = \sum_{j} w'_{ij}\big[ h_i - h_j \big] w_{ij}' = \beta S_{ij} \bar{G}, ~\frac{dw_{ij}'}{dz} = \frac{w_{ij}'}{\Delta Z}=\beta g_{ij} \bar{G}. \end{equation} -where $\beta$ is the turbulent mixing parameter or thermal transfer coefficient and $\bar{G}$ is the average mass flux of the adjacent subchannels. The $\beta$ term is the tuning parameter for the mixing model. Physically, it is a non-dimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. It is used to model the effect of the unresolved scales of motion that are produced through the averaging process. In single-phase flow no net mass exchange occurs, both momentum and energy are exchanged between subchannels, and their rates of exchange are characterized in terms of hypothetical turbulent interchange flow rates ($w_{ij}^{'H},w_{ij}^{'M}$) [!cite](TODREAS), for enthalpy and momentum respectively. The approximation that the rate of turbulent exchange for energy and momentum are equal is adopted: $w'_{ij} = w_{ij}^{'H} = w_{ij}^{'M}$. +where $\beta$ is the turbulent mixing parameter or thermal transfer coefficient and $\bar{G}$ is the average mass flux of the adjacent subchannels. The $\beta$ term is the tuning parameter for the mixing model. Physically, it is a non-dimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. It is used to model the effect of the unresolved scales of motion that are produced through the averaging process. In single-phase flow no net mass exchange occurs, both momentum and energy are exchanged between subchannels, and their rates of exchange are characterized in terms of hypothetical turbulent interchange flow rates ($w_{ij}^{'H},w_{ij}^{'M}$) [!cite](TODREAS), for enthalpy and momentum respectively. The approximation that the rate of turbulent exchange for energy and momentum are related as follows is adopted: $w'_{ij} = w_{ij}^{'H} = w_{ij}^{'M} / C_T$. Additional turbulent mixing parameters are implemented as follows: @@ -136,7 +136,6 @@ After calibrating the turbulent diffusion coefficient $\beta$ we turned our atte For quadrilateral assemblies: $C_{T} = 2.6$, $\beta = 0.006$ [!cite](kyriakopoulos2022development). - ## Discretization The collocated discretization of the variables is presented in [fig:dis] . $i,j$ are the subchannel indexes. $ij$ is the name of the gap between subchannels $i,j$. $k$ is the index in the axial direction. @@ -269,11 +268,11 @@ This is the default algorithm, where the unknown flow variables are calculated i #### Implicit segregated -In this case, the governing equations are recast in matrix form and the flow variables are calculated by solving the corresponding system. This means that variables are retrieved concurrently for the whole block. Otherwise, the solution algorithm is the same as in the default method. +In this case, the governing mass, axial momentum and crossflow momentum, equations are recast in matrix form and the flow variables are calculated by solving the corresponding system. This means that variables are retrieved concurrently for the whole block. Otherwise, the solution algorithm is the same as in the default explicit method. #### Implicit -In this case, the conservation equations are recast in matrix form and combined into a single system. The user can decide whether or not they will include the enthalpy calculation in the matrix formulation. The flow variables are calculated by solving that big system to retrieve all the unknows at the same time instead of one by one, and on all the nodes of the block: $\vec{\dot{m}}, \vec{P}, \vec{w_{ij}}, \vec{h}$. The solution algorithm is the same as in the default method, but the solver used in this version is a fixed point iteration instead of a Newton method. The system looks like this: +In this case, the governing mass, axial momentum and crossflow momentum conservation equations are recast in matrix form and combined into a single system. The system of all the subchannel equations looks like this: \begin{equation} \begin{bmatrix} @@ -297,5 +296,93 @@ In this case, the conservation equations are recast in matrix form and combined \end{bmatrix} \end{equation} -As soon as the big matrix is constructed, the solver will calculate cross-flow resistances to achieve realizable solutions. This is done by an initial calculation of the axial mass flows, crossflows, and the sum of crossflows. -The defined relaxation factor is the average crossflow divided by the maximum and added to 0.5. The added cross-flow resistance is a blending of a current value and the previous one using the relaxation factor calculated above. The current value is the maximum of the sum of cross-flows per subchannel over the minimum axial mass-flow rate. The added cross-flow resistance is added to the diagonal of $M_{ww}$ matrix. \ No newline at end of file +Since the enthalpy governing equations are uncoupled from the other equations in this otherwise monolithic system (enthalpy is coupled to the flow equations via the fluid properties update), it makes sense to lag the enthalpy solution and solve for it separately. The flow variables are calculated by solving that big system (without the enthalpy) to retrieve all the unknowns at the same time instead of one by one, and on all the nodes of the block: $\vec{\dot{m}}, \vec{P}, \vec{w_{ij}}$. The solution algorithm is the same as in the default method and the solver used is PETSc KSPSolve. + +As soon as the big matrix is constructed, the solver will calculate cross-flow resistances to maintain realizability. A distinctive feature of this method is the introduction of a *weak relaxation* logic that stabilizes and accelerates convergence of the coupled $mass flow: (\dot{\mathbf{m}})$, $pressure: (\mathbf{P})$, and $crossflow:(\mathbf{w}_{ij})$ fields in a $Q{=}3$ block-nested linear system with matrix blocks $M_{ij}$ and right-hand-side blocks $\mathbf{b}_i$ that represent the individual governing equations. Note that the solution is influenced by the stabilization method and its coefficients. + +#### 1. Fast scale estimates + +From the axial- and cross-momentum rows, the code forms quick, diagonally preconditioned estimates: +\begin{equation} +\begin{aligned} +\hat{\mathbf m} &= M_{pm}\,\mathbf m, \\ +\hat{\mathbf p} &= \frac{\hat{\mathbf m}}{\operatorname{diag}(M_{pp}) + \varepsilon_p\mathbf 1},\\ +\hat{\mathbf W} &= \frac{M_{wp}\,\hat{\mathbf p} - \mathbf b_w}{\operatorname{diag}(M_{ww}) + \varepsilon_W\mathbf 1}, +\end{aligned} +\end{equation} +with small safeguards $\varepsilon_p,\varepsilon_W\sim 10^{-10}$ to avoid division by zero. Using $\hat{\mathbf W}$, the per-channel crossflow sum $\sum_{j} w_{ij}$ is assembled into a vector $\mathrm{sumw_{ij}}_{\mathrm{loc}}$. + +#### 2. Crossflow relaxation parameter + +Two guarded scalars are computed: + +\begin{equation} +\begin{aligned} +m_{\min} &= \max\big(\min |\mathbf m|,\; 10^{-10}\big),\\ +S_{\max} &= \max\Big(\max |\mathrm{sumw_{ijloc}}|,\; 10^{-10}\Big) +\end{aligned} +\end{equation} + +Additionally, a mean inter-iteration change for crossflow is formed +\begin{equation} +r_{\mathrm{base}} = \operatorname{mean}\big(\big|\mathbf W^{(k)}| - |\mathbf W^{(k-1)}\big|\big), +\end{equation} +leading to a relaxation factor +\begin{equation} +r = \frac{r_{\mathrm{base}}}{\max(S_{\max}, \varepsilon)} + 0.5,\qquad \varepsilon\sim10^{-10}. +\end{equation} +The +0.5 offset biases toward mild under-relaxation. + +#### 3. Crossflow resistance inflation + +A cross-coupling resistance is estimated and smoothed: +\begin{equation} +\begin{aligned} +\tilde K &= \frac{S_{\max}}{m_{\min}}, & +K^\star &= 0.9\,\tilde K + 0.1\,K_{\text{old}}, & +K &= r\,K^\star. +\end{aligned} +\end{equation} +After smoothing, the provisional crossflow resistance $K$ is mapped through a piecewise lower-bound function that enforces minimum safe damping levels in specific ranges. + +\begin{equation} +K \rightarrow +\begin{cases} +K , & K >= 10, \\ +1.0, & 1 \leq K < 10, \\ +0.5, & 0.1 \leq K < 1, \\ +\frac{1}{3}, & 0.01 \leq K < 0.1, \\ +0.1, & 0.001 \leq K < 0.01, \\ +K, & K < 10^{-3}. +\end{cases} +\end{equation} + +This mapping acts as a {snap-up} rule for the crossflow resistance $K$ over the range $[10^{-3}, 10]$: +it raises $K$ out of weak-damping intervals but leaves very small and very large +values unchanged. The purpose is to maintain numerical stability and adequate +diagonal dominance in the cross-momentum equations without introducing full quantization or "bucketing". + +Finally, $K$ is added to the diagonal of the cross-momentum block, +\begin{equation} +M_{ww} \;\leftarrow\; M_{ww} + K\,I, +\end{equation} +thereby increasing diagonal dominance and improving conditioning for the crossflow equations. Note that this treatment does influence the cross-flow distribution solution. + +#### 4. Per-equation under-relaxation + +Classical linear under-relaxation is applied automatically and separately to each equation $f\in\{\mathbf m,\mathbf p,\mathbf W\}$ using factors +\begin{equation} +\alpha_m=1.0,\qquad \alpha_p=1.0,\qquad \alpha_W=0.1. +\end{equation} +For each equation, with the corresponding diagonal $D_f=\operatorname{diag}(M_{ff})$, the system is modified as +\begin{equation} +\begin{aligned} +M_{ff} &\leftarrow \frac{1}{\alpha_f} M_{ff}, \\ +\mathbf b_f &\leftarrow \mathbf b_f + (1-\alpha_f)\,D_f\,\mathbf x^{\text{old}}_f. +\end{aligned} +\end{equation} +This standard construction ensures that solving the modified linear system yields the *under-relaxed* update for equation $f$. In practice, only $\mathbf W$ is strongly damped, while $\mathbf m$ and $\mathbf p$ can be solved without additional damping. This relaxation happens inside the temperature loop. + +#### 5. Net effect + +The combination of (i) safeguarded scale estimation, (ii) adaptive, time smoothed, and piecewise snapped added crossflow resistance, and (iii) selective under-relaxation produces a more diagonally dominant and robust nested solve that tolerates rapid changes in crossflow while preserving good convergence properties for mass flow and pressure. diff --git a/modules/subchannel/examples/MultiApp/fuel_assembly.i b/modules/subchannel/examples/MultiApp/fuel_assembly.i index 5440371e0015..1ae88067bef9 100644 --- a/modules/subchannel/examples/MultiApp/fuel_assembly.i +++ b/modules/subchannel/examples/MultiApp/fuel_assembly.i @@ -145,7 +145,6 @@ duct_inside = '${fparse duct_outside - 2 * duct_thickness}' implicit = false segregated = true staggered_pressure = false - monolithic_thermal = false # Tolerances P_tol = 1.0e-4 diff --git a/modules/subchannel/examples/duct/test.i b/modules/subchannel/examples/duct/test.i index ae88e24723ed..7c88bacf7360 100644 --- a/modules/subchannel/examples/duct/test.i +++ b/modules/subchannel/examples/duct/test.i @@ -92,7 +92,6 @@ P_out = 2.0e5 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_multiapps = true verbose_subchannel = false [] diff --git a/modules/subchannel/include/problems/SubChannel1PhaseProblem.h b/modules/subchannel/include/problems/SubChannel1PhaseProblem.h index e33d8d200c02..0127d45ca6dd 100644 --- a/modules/subchannel/include/problems/SubChannel1PhaseProblem.h +++ b/modules/subchannel/include/problems/SubChannel1PhaseProblem.h @@ -138,6 +138,22 @@ class SubChannel1PhaseProblem : public ExternalProblem, public PostprocessorInte } } + /** + * Solve a linear system (A * x = rhs) with a simple PCJACOBI KSP and populate the + * enthalpy solution into _h_soln for nodes [first_node, last_node]. + * + * Uses member tolerances (_rtol, _atol, _dtol, _maxit), mesh (_subchannel_mesh), + * channel count (_n_channels), and error/solution handles (mooseError, _h_soln). + * + * @param A PETSc matrix (operators) + * @param rhs PETSc vector (right-hand side) + * @param first_node inclusive start axial node index + * @param last_node inclusive end axial node index + * @param ksp_prefix options prefix for KSP (e.g. "h_sys_"), may be nullptr + */ + PetscErrorCode solveAndPopulateEnthalpy( + Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix); + PetscErrorCode cleanUp(); SubChannelMesh & _subchannel_mesh; /// number of axial blocks @@ -201,8 +217,6 @@ class SubChannel1PhaseProblem : public ExternalProblem, public PostprocessorInte const bool _staggered_pressure_bool; /// Segregated solve const bool _segregated_bool; - /// Thermal monolithic bool - const bool _monolithic_thermal_bool; /// Boolean to printout information related to subchannel solve const bool _verbose_subchannel; /// Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement, Dpin diff --git a/modules/subchannel/src/problems/QuadSubChannel1PhaseProblem.C b/modules/subchannel/src/problems/QuadSubChannel1PhaseProblem.C index 53a96ddc6b97..130d61c65e09 100644 --- a/modules/subchannel/src/problems/QuadSubChannel1PhaseProblem.C +++ b/modules/subchannel/src/problems/QuadSubChannel1PhaseProblem.C @@ -806,43 +806,8 @@ QuadSubChannel1PhaseProblem::computeh(int iblock) LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_cross_derivative_rhs)); LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_added_heat_rhs)); - if (_segregated_bool || (!_monolithic_thermal_bool)) - { - // Assembly the matrix system - KSP ksploc; - PC pc; - Vec sol; - LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol)); - LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc)); - LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat)); - LibmeshPetscCall(KSPGetPC(ksploc, &pc)); - LibmeshPetscCall(PCSetType(pc, PCJACOBI)); - LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit)); - LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_")); - LibmeshPetscCall(KSPSetFromOptions(ksploc)); - LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol)); - PetscScalar * xx; - LibmeshPetscCall(VecGetArray(sol, &xx)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) - { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); - auto h_out = xx[iz_ind * _n_channels + i_ch]; - if (h_out < 0) - { - mooseError(name(), - " : Calculation of negative Enthalpy h_out = : ", - h_out, - " Axial Level= : ", - iz); - } - _h_soln->set(node_out, h_out); - } - } - LibmeshPetscCall(KSPDestroy(&ksploc)); - LibmeshPetscCall(VecDestroy(&sol)); - } + // Use system to solve for and populate enthalpy + LibmeshPetscCall(this->solveAndPopulateEnthalpy( + _hc_sys_h_mat, _hc_sys_h_rhs, first_node, last_node, "h_sys_")); } } diff --git a/modules/subchannel/src/problems/SubChannel1PhaseProblem.C b/modules/subchannel/src/problems/SubChannel1PhaseProblem.C index c0d8060f5a78..0d594a1507be 100644 --- a/modules/subchannel/src/problems/SubChannel1PhaseProblem.C +++ b/modules/subchannel/src/problems/SubChannel1PhaseProblem.C @@ -90,8 +90,6 @@ SubChannel1PhaseProblem::validParams() "staggered_pressure", false, "Boolean to define the use of explicit or implicit solution."); params.addParam( "segregated", true, "Boolean to define whether to use a segregated solution."); - params.addParam( - "monolithic_thermal", false, "Boolean to define whether to use thermal monolithic solve."); params.addParam( "verbose_subchannel", false, "Boolean to print out information related to subchannel solve."); params.addParam( @@ -148,7 +146,6 @@ SubChannel1PhaseProblem::SubChannel1PhaseProblem(const InputParameters & params) _implicit_bool(getParam("implicit")), _staggered_pressure_bool(getParam("staggered_pressure")), _segregated_bool(getParam("segregated")), - _monolithic_thermal_bool(getParam("monolithic_thermal")), _verbose_subchannel(getParam("verbose_subchannel")), _deformation(getParam("deformation")), _fp(nullptr), @@ -2011,6 +2008,56 @@ SubChannel1PhaseProblem::petscSnesSolver(int iblock, PetscFunctionReturn(LIBMESH_PETSC_SUCCESS); } +PetscErrorCode +SubChannel1PhaseProblem::solveAndPopulateEnthalpy( + Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix) +{ + PetscFunctionBegin; + + // Create solution vector with rhs layout + Vec x = nullptr; + LibmeshPetscCall(VecDuplicate(rhs, &x)); + + // KSP setup + KSP ksp = nullptr; + PC pc = nullptr; + LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); + LibmeshPetscCall(KSPSetOperators(ksp, A, A)); + LibmeshPetscCall(KSPGetPC(ksp, &pc)); + LibmeshPetscCall(PCSetType(pc, PCJACOBI)); + LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit)); + if (ksp_prefix && *ksp_prefix) + LibmeshPetscCall(KSPSetOptionsPrefix(ksp, ksp_prefix)); + LibmeshPetscCall(KSPSetFromOptions(ksp)); + + // Solve + LibmeshPetscCall(KSPSolve(ksp, rhs, x)); + + // Scatter to _h_soln with sanity checks + PetscScalar * xx = nullptr; + LibmeshPetscCall(VecGetArray(x, &xx)); + for (unsigned int iz = first_node; iz <= last_node; ++iz) + { + const unsigned int iz_ind = iz - first_node; + for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch) + { + auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); + const PetscScalar h_out = xx[iz_ind * _n_channels + i_ch]; + if (h_out < 0.0) + mooseError( + name(), " : Calculation of negative Enthalpy h_out = ", h_out, " Axial Level = ", iz); + _h_soln->set(node_out, h_out); + } + } + LibmeshPetscCall(VecRestoreArray(x, &xx)); + + // Cleanup + LibmeshPetscCall(KSPDestroy(&ksp)); + LibmeshPetscCall(VecDestroy(&x)); + + PetscFunctionReturn(LIBMESH_PETSC_SUCCESS); +} + Real SubChannel1PhaseProblem::computeAddedHeatDuct(unsigned int i_ch, unsigned int iz) { @@ -2044,217 +2091,144 @@ SubChannel1PhaseProblem::computeAddedHeatDuct(unsigned int i_ch, unsigned int iz PetscErrorCode SubChannel1PhaseProblem::implicitPetscSolve(int iblock) { - bool lag_block_thermal_solve = true; - Vec b_nest, x_nest; /* approx solution, RHS, exact solution */ - Mat A_nest; /* linear system matrix */ - KSP ksp; /* linear solver context */ - PC pc; /* preconditioner context */ - PetscFunctionBegin; - PetscInt Q = _monolithic_thermal_bool ? 4 : 3; - std::vector mat_array(Q * Q); - std::vector vec_array(Q); - - /// Initializing flags - bool _axial_mass_flow_tight_coupling = true; - bool _pressure_axial_momentum_tight_coupling = true; - bool _pressure_cross_momentum_tight_coupling = true; - unsigned int first_node = iblock * _block_size + 1; - unsigned int last_node = (iblock + 1) * _block_size; - - /// Assembling matrices - // Computing sum of crossflows with previous iteration - computeSumWij(iblock); - // Assembling axial flux matrix - computeMdot(iblock); - // Computing turbulent crossflow with previous step axial mass flows - computeWijPrime(iblock); - // Assembling for Pressure Drop matrix - computeDP(iblock); - // Assembling pressure matrix - computeP(iblock); - // Assembling cross fluxes matrix - computeWijResidual(iblock); - // If monolithic solve - Assembling enthalpy matrix - if (_monolithic_thermal_bool) - computeh(iblock); - - if (_verbose_subchannel) + // ---------- helper functions ----------------------------- + auto V = [&](const std::string & s) { - _console << "Starting nested system." << std::endl; - _console << "Number of simultaneous variables: " << Q << std::endl; - } - // Mass conservation - PetscInt field_num = 0; - LibmeshPetscCall( - MatDuplicate(_mc_axial_convection_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 0])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY)); - mat_array[Q * field_num + 1] = NULL; - if (_axial_mass_flow_tight_coupling) - { - LibmeshPetscCall(MatDuplicate(_mc_sumWij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY)); - } - else - { - mat_array[Q * field_num + 2] = NULL; - } - if (_monolithic_thermal_bool) - { - mat_array[Q * field_num + 3] = NULL; - } - LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &vec_array[field_num])); - LibmeshPetscCall(VecCopy(_mc_axial_convection_rhs, vec_array[field_num])); - if (!_axial_mass_flow_tight_coupling) + if (_verbose_subchannel) + _console << s << std::endl; + }; + + auto DupMatAssembled = [&](Mat src, Mat * dst) { - Vec sumWij_loc; - LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sumWij_loc)); - LibmeshPetscCall(VecSet(sumWij_loc, 0.0)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) + if (src) { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); - - PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out); - PetscInt row_vec_2 = i_ch + _n_channels * iz_ind; - LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec_2, &value_vec_2, ADD_VALUES)); - } + LibmeshPetscCall(MatDuplicate(src, MAT_COPY_VALUES, dst)); + LibmeshPetscCall(MatAssemblyBegin(*dst, MAT_FINAL_ASSEMBLY)); + LibmeshPetscCall(MatAssemblyEnd(*dst, MAT_FINAL_ASSEMBLY)); } - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sumWij_loc)); - LibmeshPetscCall(VecDestroy(&sumWij_loc)); - } + else + *dst = NULL; + }; - if (_verbose_subchannel) - _console << "Mass ok." << std::endl; - // Axial momentum conservation - field_num = 1; - if (_pressure_axial_momentum_tight_coupling) - { - LibmeshPetscCall( - MatDuplicate(_amc_sys_mdot_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 0])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY)); - } - else - { - mat_array[Q * field_num + 0] = NULL; - } - LibmeshPetscCall( - MatDuplicate(_amc_pressure_force_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 1])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY)); - mat_array[Q * field_num + 2] = NULL; - if (_monolithic_thermal_bool) + auto DupVecCopy = [&](Vec src, Vec * dst) { - mat_array[Q * field_num + 3] = NULL; - } - LibmeshPetscCall(VecDuplicate(_amc_pressure_force_rhs, &vec_array[field_num])); - LibmeshPetscCall(VecCopy(_amc_pressure_force_rhs, vec_array[field_num])); - if (_pressure_axial_momentum_tight_coupling) - { - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _amc_sys_mdot_rhs)); - } - else - { - unsigned int last_node = (iblock + 1) * _block_size; - unsigned int first_node = iblock * _block_size + 1; - LibmeshPetscCall(populateVectorFromHandle( - _prod, *_mdot_soln, first_node, last_node, _n_channels)); - Vec ls; - LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &ls)); - LibmeshPetscCall(MatMult(_amc_sys_mdot_mat, _prod, ls)); - LibmeshPetscCall(VecAXPY(ls, -1.0, _amc_sys_mdot_rhs)); - LibmeshPetscCall(VecAXPY(vec_array[field_num], -1.0, ls)); - LibmeshPetscCall(VecDestroy(&ls)); - } + LibmeshPetscCall(VecDuplicate(src, dst)); + LibmeshPetscCall(VecCopy(src, *dst)); + }; - if (_verbose_subchannel) - _console << "Lin mom OK." << std::endl; + const PetscInt Q = 3; // [mass conservation, axial momentum, cross momentum] - // Cross momentum conservation - field_num = 2; - mat_array[Q * field_num + 0] = NULL; - if (_pressure_cross_momentum_tight_coupling) - { - LibmeshPetscCall( - MatDuplicate(_cmc_pressure_force_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 1])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY)); - } - else - { - mat_array[Q * field_num + 1] = NULL; - } + // small indexer + auto Idx = [&](PetscInt r, PetscInt c) { return Q * r + c; }; - LibmeshPetscCall(MatDuplicate(_cmc_sys_Wij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2])); - LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY)); - LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY)); - if (_monolithic_thermal_bool) - { - mat_array[Q * field_num + 3] = NULL; - } + // arrays that MUST be declared before lambdas use them + std::vector mat_array(Q * Q, NULL); + std::vector vec_array(Q, NULL); - LibmeshPetscCall(VecDuplicate(_cmc_sys_Wij_rhs, &vec_array[field_num])); - LibmeshPetscCall(VecCopy(_cmc_sys_Wij_rhs, vec_array[field_num])); - if (_pressure_cross_momentum_tight_coupling) + // generic assembler for one governing equation row in the nested matrix + auto AssembleEquation = [&](PetscInt f, + Mat A0, + Mat A1, + Mat A2, // three blocks in row f (can be nullptr) + Vec rhs, // base RHS for equation f + Vec rhs_add, // optional extra RHS to add (can be nullptr) + const char * label) // e.g. "Mass", "Lin mom", "Cross mom" { - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _cmc_pressure_force_rhs)); - } - else + DupMatAssembled(A0, &mat_array[Idx(f, 0)]); + DupMatAssembled(A1, &mat_array[Idx(f, 1)]); + DupMatAssembled(A2, &mat_array[Idx(f, 2)]); + DupVecCopy(rhs, &vec_array[f]); + if (rhs_add) + LibmeshPetscCall(VecAXPY(vec_array[f], 1.0, rhs_add)); + V(std::string(label) + " system assembled"); + }; + + // ----------------------------------------------------------------------------- + // Helper lambda that applies per-equation under-relaxation by modifying BOTH the + // matrix block and the RHS for that equation. + // + // Specifically, with A_ff for the equation and D = diag(A_ff): + // 1) Matrix diagonal scaling: D <- D / alpha, then A_ff's diagonal is replaced + // with D. For alpha < 1 this increases diagonal dominance (more damping). + // 2) RHS blending with the previous solution x_old: + // rhs_f <- rhs_f + (1 - alpha) * (D / alpha) * x_old + // where x_old is provided by the caller via `populate(work)`. + // + // Net effect: the solved x satisfies + // A_ff x = rhs_f_original + ((1 - alpha)/alpha) * D * (x_old - x), + // which damps updates toward x_old without changing the converged solution. + // ----------------------------------------------------------------------------- + auto RelaxEquation = + [&](Mat A_ff, Vec rhs_f, Vec like_vec, Vec work, PetscScalar alpha, auto && populate) { - Vec sol_holder_P; - LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps)); - LibmeshPetscCall(populateVectorFromHandle( - _prodp, *_P_soln, iblock * _block_size, (iblock + 1) * _block_size - 1, _n_channels)); + Vec d = nullptr; + LibmeshPetscCall(VecDuplicate(like_vec, &d)); - LibmeshPetscCall(MatMult(_cmc_pressure_force_mat, _prodp, sol_holder_P)); - LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs)); - LibmeshPetscCall(VecScale(sol_holder_P, 1.0)); - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sol_holder_P)); - LibmeshPetscCall(VecDestroy(&sol_holder_P)); - } + // 1) A_ff: diag <- diag / alpha + LibmeshPetscCall(MatGetDiagonal(A_ff, d)); + LibmeshPetscCall(VecScale(d, 1.0 / alpha)); + LibmeshPetscCall(MatDiagonalSet(A_ff, d, INSERT_VALUES)); - if (_verbose_subchannel) - _console << "Cross mom ok." << std::endl; + // 2) work <- x_old (caller-provided populator) + LibmeshPetscCall(populate(work)); - // 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])); - if (lag_block_thermal_solve) - { - LibmeshPetscCall(VecZeroEntries(vec_array[field_num])); - LibmeshPetscCall(VecShift(vec_array[field_num], 1.0)); - } + // 3) rhs_f += (1 - alpha) * (diag .* work) + LibmeshPetscCall(VecScale(d, (1.0 - alpha))); + LibmeshPetscCall(VecPointwiseMult(work, work, d)); + LibmeshPetscCall(VecAXPY(rhs_f, 1.0, work)); - if (_verbose_subchannel) - _console << "Energy ok." << std::endl; - } + LibmeshPetscCall(VecDestroy(&d)); + }; + + // indices + const unsigned int first_node = iblock * _block_size + 1; + const unsigned int last_node = (iblock + 1) * _block_size; + + // ---------- assemble per-block operators ----------------- + computeSumWij(iblock); + computeMdot(iblock); + computeWijPrime(iblock); + computeDP(iblock); + computeP(iblock); + computeWijResidual(iblock); - // Relaxing linear system - // Weaker relaxation + V("Starting nested system."); + + // Populate nested matrix with the individual physics + // equation 0: Mass conservation + AssembleEquation(/*f=*/0, + /*A0=*/_mc_axial_convection_mat, + /*A1=*/nullptr, + /*A2=*/_mc_sumWij_mat, + /*rhs=*/_mc_axial_convection_rhs, + /*rhs_add=*/nullptr, + /*label=*/"Mass"); + + // equation 1: Axial momentum conservation + AssembleEquation(/*f=*/1, + /*A0=*/_amc_sys_mdot_mat, + /*A1=*/_amc_pressure_force_mat, + /*A2=*/nullptr, + /*rhs=*/_amc_pressure_force_rhs, + /*rhs_add=*/_amc_sys_mdot_rhs, + /*label=*/"Lin mom"); + + // equation 2: Cross momentum conservation + AssembleEquation(/*f=*/2, + /*A0=*/nullptr, + /*A1=*/_cmc_pressure_force_mat, + /*A2=*/_cmc_sys_Wij_mat, + /*rhs=*/_cmc_sys_Wij_rhs, + /*rhs_add=*/_cmc_pressure_force_rhs, + /*label=*/"Cross mom"); + + // ========================== Relaxation ==================== if (true) { - // Estimating cross-flow resistances to achieve realizable solves LibmeshPetscCall(populateVectorFromHandle( _prod, *_mdot_soln, first_node, last_node, _n_channels)); + Vec mdot_estimate; LibmeshPetscCall(createPetscVector(mdot_estimate, _block_size * _n_channels)); Vec pmat_diag; @@ -2266,10 +2240,6 @@ SubChannel1PhaseProblem::implicitPetscSolve(int iblock) LibmeshPetscCall(VecSet(unity_vec, 1.0)); Vec sol_holder_P; LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps)); - Vec diag_Wij_loc; - LibmeshPetscCall(createPetscVector(diag_Wij_loc, _block_size * _n_gaps)); - Vec Wij_estimate; - LibmeshPetscCall(createPetscVector(Wij_estimate, _block_size * _n_gaps)); Vec unity_vec_Wij; LibmeshPetscCall(createPetscVector(unity_vec_Wij, _block_size * _n_gaps)); LibmeshPetscCall(VecSet(unity_vec_Wij, 1.0)); @@ -2277,21 +2247,27 @@ SubChannel1PhaseProblem::implicitPetscSolve(int iblock) LibmeshPetscCall(createPetscVector(_Wij_loc_vec, _block_size * _n_gaps)); Vec _Wij_old_loc_vec; LibmeshPetscCall(createPetscVector(_Wij_old_loc_vec, _block_size * _n_gaps)); - LibmeshPetscCall(MatMult(mat_array[Q], _prod, mdot_estimate)); + + // ---- scale estimates ---- + // mdot_estimate = A(1,0) * mdot + LibmeshPetscCall(MatMult(mat_array[Q /* (1,0) */], _prod, mdot_estimate)); + + // p_estimate = mdot_est / (diag(A(1,1)) + eps) LibmeshPetscCall(MatGetDiagonal(mat_array[Q + 1], pmat_diag)); LibmeshPetscCall(VecAXPY(pmat_diag, 1e-10, unity_vec)); LibmeshPetscCall(VecPointwiseDivide(p_estimate, mdot_estimate, pmat_diag)); + + // sol_holder_P = A(2,1) * p_estimate - rhs_cmc_pressure LibmeshPetscCall(MatMult(mat_array[2 * Q + 1], p_estimate, sol_holder_P)); LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs)); - LibmeshPetscCall(MatGetDiagonal(mat_array[2 * Q + 2], diag_Wij_loc)); - LibmeshPetscCall(VecAXPY(diag_Wij_loc, 1e-10, unity_vec_Wij)); - LibmeshPetscCall(VecPointwiseDivide(Wij_estimate, sol_holder_P, diag_Wij_loc)); + + // sumWij_loc from sol_holder_P (accumulate) Vec sumWij_loc; LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) + for (unsigned int iz = first_node; iz <= last_node; ++iz) { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) + const auto iz_ind = iz - first_node; + for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch) { PetscScalar sumWij = 0.0; unsigned int counter = 0; @@ -2309,18 +2285,19 @@ SubChannel1PhaseProblem::implicitPetscSolve(int iblock) LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec, &sumWij, INSERT_VALUES)); } } + LibmeshPetscCall(VecAssemblyBegin(sumWij_loc)); + LibmeshPetscCall(VecAssemblyEnd(sumWij_loc)); + // ---- robust scale measurements ---- PetscScalar min_mdot; LibmeshPetscCall(VecAbs(_prod)); LibmeshPetscCall(VecMin(_prod, NULL, &min_mdot)); - if (_verbose_subchannel) - _console << "Minimum estimated mdot: " << min_mdot << std::endl; + V("Minimum estimated mdot: " + std::to_string(min_mdot)); LibmeshPetscCall(VecAbs(sumWij_loc)); LibmeshPetscCall(VecMax(sumWij_loc, NULL, &_max_sumWij)); _max_sumWij = std::max(1e-10, _max_sumWij); - if (_verbose_subchannel) - _console << "Maximum estimated Wij: " << _max_sumWij << std::endl; + V("Maximum estimated Wij: " + std::to_string(_max_sumWij)); LibmeshPetscCall(populateVectorFromDense>( _Wij_loc_vec, _Wij, first_node, last_node, _n_gaps)); @@ -2329,6 +2306,7 @@ SubChannel1PhaseProblem::implicitPetscSolve(int iblock) _Wij_old_loc_vec, _Wij_old, first_node, last_node, _n_gaps)); LibmeshPetscCall(VecAbs(_Wij_old_loc_vec)); LibmeshPetscCall(VecAXPY(_Wij_loc_vec, -1.0, _Wij_old_loc_vec)); + PetscScalar relax_factor; LibmeshPetscCall(VecAbs(_Wij_loc_vec)); #if !PETSC_VERSION_LESS_THAN(3, 16, 0) @@ -2338,305 +2316,192 @@ SubChannel1PhaseProblem::implicitPetscSolve(int iblock) relax_factor /= _block_size * _n_gaps; #endif relax_factor = relax_factor / _max_sumWij + 0.5; - if (_verbose_subchannel) - _console << "Relax base value: " << relax_factor << std::endl; + V("Relax base value: " + std::to_string(relax_factor)); - PetscScalar resistance_relaxation = 0.9; + // ---- crossflow resistance inflation ---- + const PetscScalar resistance_relaxation = 0.9; _added_K = _max_sumWij / min_mdot; - if (_verbose_subchannel) - _console << "New cross resistance: " << _added_K << std::endl; + V("New cross resistance: " + std::to_string(_added_K)); _added_K = (_added_K * resistance_relaxation + (1.0 - resistance_relaxation) * _added_K_old) * relax_factor; - if (_verbose_subchannel) - _console << "Relaxed cross resistance: " << _added_K << std::endl; + V("Relaxed cross resistance: " + std::to_string(_added_K)); + + // Snap-up lower-bounding if (_added_K < 10 && _added_K >= 1.0) - _added_K = 1.0; //(1.0 - resistance_relaxation); + _added_K = 1.0; if (_added_K < 1.0 && _added_K >= 0.1) _added_K = 0.5; if (_added_K < 0.1 && _added_K >= 0.01) _added_K = 1. / 3.; if (_added_K < 1e-2 && _added_K >= 1e-3) _added_K = 0.1; - if (_added_K < 1e-3) - _added_K = 1.0 * _added_K; - if (_verbose_subchannel) - _console << "Actual added cross resistance: " << _added_K << std::endl; + V("Actual added cross resistance: " + std::to_string(_added_K)); LibmeshPetscCall(VecScale(unity_vec_Wij, _added_K)); _added_K_old = _added_K; - // Adding cross resistances LibmeshPetscCall(MatDiagonalSet(mat_array[2 * Q + 2], unity_vec_Wij, ADD_VALUES)); + + // ---- cleanup temp vectors used above ---- LibmeshPetscCall(VecDestroy(&mdot_estimate)); LibmeshPetscCall(VecDestroy(&pmat_diag)); LibmeshPetscCall(VecDestroy(&unity_vec)); LibmeshPetscCall(VecDestroy(&p_estimate)); LibmeshPetscCall(VecDestroy(&sol_holder_P)); - LibmeshPetscCall(VecDestroy(&diag_Wij_loc)); LibmeshPetscCall(VecDestroy(&unity_vec_Wij)); - LibmeshPetscCall(VecDestroy(&Wij_estimate)); LibmeshPetscCall(VecDestroy(&sumWij_loc)); LibmeshPetscCall(VecDestroy(&_Wij_loc_vec)); LibmeshPetscCall(VecDestroy(&_Wij_old_loc_vec)); - // Auto-computing relaxation factors - PetscScalar relaxation_factor_mdot, relaxation_factor_P, relaxation_factor_Wij; - relaxation_factor_mdot = 1.0; - relaxation_factor_P = 1.0; // std::exp(-5.0); - relaxation_factor_Wij = 0.1; - - if (_verbose_subchannel) - { - _console << "Relax mdot: " << relaxation_factor_mdot << std::endl; - _console << "Relax P: " << relaxation_factor_P << std::endl; - _console << "Relax Wij: " << relaxation_factor_Wij << std::endl; - } - - PetscInt field_num = 0; - Vec diag_mdot; - LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_mdot)); - LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_mdot)); - LibmeshPetscCall(VecScale(diag_mdot, 1.0 / relaxation_factor_mdot)); - LibmeshPetscCall( - MatDiagonalSet(mat_array[Q * field_num + field_num], diag_mdot, INSERT_VALUES)); - LibmeshPetscCall(populateVectorFromHandle( - _prod, *_mdot_soln, first_node, last_node, _n_channels)); - LibmeshPetscCall(VecScale(diag_mdot, (1.0 - relaxation_factor_mdot))); - LibmeshPetscCall(VecPointwiseMult(_prod, _prod, diag_mdot)); - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _prod)); - LibmeshPetscCall(VecDestroy(&diag_mdot)); - - if (_verbose_subchannel) - _console << "mdot relaxed" << std::endl; - - field_num = 1; - Vec diag_P; - LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_P)); - LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_P)); - LibmeshPetscCall(VecScale(diag_P, 1.0 / relaxation_factor_P)); - LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_P, INSERT_VALUES)); - if (_verbose_subchannel) - _console << "Mat assembled" << std::endl; - LibmeshPetscCall(populateVectorFromHandle( - _prod, *_P_soln, first_node, last_node, _n_channels)); - LibmeshPetscCall(VecScale(diag_P, (1.0 - relaxation_factor_P))); - LibmeshPetscCall(VecPointwiseMult(_prod, _prod, diag_P)); - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _prod)); - LibmeshPetscCall(VecDestroy(&diag_P)); - - if (_verbose_subchannel) - _console << "P relaxed" << std::endl; - - field_num = 2; - Vec diag_Wij; - LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_Wij)); - LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_Wij)); - LibmeshPetscCall(VecScale(diag_Wij, 1.0 / relaxation_factor_Wij)); - LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_Wij, INSERT_VALUES)); - LibmeshPetscCall(populateVectorFromDense>( - _Wij_vec, _Wij, first_node, last_node, _n_gaps)); - LibmeshPetscCall(VecScale(diag_Wij, (1.0 - relaxation_factor_Wij))); - LibmeshPetscCall(VecPointwiseMult(_Wij_vec, _Wij_vec, diag_Wij)); - LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _Wij_vec)); - LibmeshPetscCall(VecDestroy(&diag_Wij)); - - if (_verbose_subchannel) - _console << "Wij relaxed" << std::endl; + // ---- per-equation under-relaxation ---- + const PetscScalar relaxation_factor_mdot = 1.0; + const PetscScalar relaxation_factor_P = 1.0; + const PetscScalar relaxation_factor_Wij = 0.1; + + V("Relax mdot: " + std::to_string(relaxation_factor_mdot)); + V("Relax P: " + std::to_string(relaxation_factor_P)); + V("Relax Wij: " + std::to_string(relaxation_factor_Wij)); + + // mdot + RelaxEquation(mat_array[Idx(0, 0)], + vec_array[0], + vec_array[0], + _prod, + relaxation_factor_mdot, + [&](Vec dst) + { + return populateVectorFromHandle( + dst, *_mdot_soln, first_node, last_node, _n_channels); + }); + V("mdot relaxed"); + + // pressure + RelaxEquation(mat_array[Idx(1, 1)], + vec_array[1], + vec_array[1], + _prod, + relaxation_factor_P, + [&](Vec dst) + { + return populateVectorFromHandle( + dst, *_P_soln, first_node, last_node, _n_channels); + }); + V("P relaxed"); + + // crossflow + RelaxEquation(mat_array[Idx(2, 2)], + vec_array[2], + vec_array[2], + _Wij_vec, + relaxation_factor_Wij, + [&](Vec dst) + { + return populateVectorFromDense>( + dst, _Wij, first_node, last_node, _n_gaps); + }); + V("Wij relaxed"); } - if (_verbose_subchannel) - _console << "Linear solver relaxed." << std::endl; + V("Linear solver relaxed"); - // Creating nested matrices + // ======================== Create and configure KSP ========================= + Mat A_nest; + Vec b_nest; + Vec x_nest; LibmeshPetscCall(MatCreateNest(PETSC_COMM_SELF, Q, NULL, Q, NULL, mat_array.data(), &A_nest)); LibmeshPetscCall(VecCreateNest(PETSC_COMM_SELF, Q, NULL, vec_array.data(), &b_nest)); - if (_verbose_subchannel) - _console << "Nested system created." << std::endl; + V("Nested system created"); - /// Setting up linear solver - // Creating linear solver + KSP ksp; + PC pc; LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksp)); LibmeshPetscCall(KSPSetType(ksp, KSPFGMRES)); - // Setting KSP operators LibmeshPetscCall(KSPSetOperators(ksp, A_nest, A_nest)); - // Set KSP and PC options LibmeshPetscCall(KSPGetPC(ksp, &pc)); LibmeshPetscCall(PCSetType(pc, PCFIELDSPLIT)); LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit)); - // Splitting fields + + // split equations std::vector rows(Q); - // IS rows[Q]; - PetscInt M = 0; LibmeshPetscCall(MatNestGetISs(A_nest, rows.data(), NULL)); for (PetscInt j = 0; j < Q; ++j) { - IS expand1; - LibmeshPetscCall(ISDuplicate(rows[M], &expand1)); - M += 1; - LibmeshPetscCall(PCFieldSplitSetIS(pc, NULL, expand1)); - LibmeshPetscCall(ISDestroy(&expand1)); + IS part; + LibmeshPetscCall(ISDuplicate(rows[j], &part)); + LibmeshPetscCall(PCFieldSplitSetIS(pc, NULL, part)); + LibmeshPetscCall(ISDestroy(&part)); } - if (_verbose_subchannel) - _console << "Linear solver assembled." << std::endl; + V("Linear solver assembled"); - /// Solving + // ============================== Solve ===================================== LibmeshPetscCall(VecDuplicate(b_nest, &x_nest)); LibmeshPetscCall(VecSet(x_nest, 0.0)); LibmeshPetscCall(KSPSolve(ksp, b_nest, x_nest)); - /// Destroying solver elements + // destroy solver containers first LibmeshPetscCall(VecDestroy(&b_nest)); LibmeshPetscCall(MatDestroy(&A_nest)); LibmeshPetscCall(KSPDestroy(&ksp)); for (PetscInt i = 0; i < Q * Q; i++) - { LibmeshPetscCall(MatDestroy(&mat_array[i])); - } for (PetscInt i = 0; i < Q; i++) - { LibmeshPetscCall(VecDestroy(&vec_array[i])); - } - if (_verbose_subchannel) - _console << "Solver elements destroyed." << std::endl; + V("Solver elements destroyed"); - /// Recovering the solutions + // ====================== Extract & scatter the solution ===================== Vec sol_mdot, sol_p, sol_Wij; - if (_verbose_subchannel) - _console << "Vectors created." << std::endl; + V("Vectors to hold solution created"); PetscInt num_vecs; Vec * loc_vecs; LibmeshPetscCall(VecNestGetSubVecs(x_nest, &num_vecs, &loc_vecs)); - if (_verbose_subchannel) - _console << "Starting extraction." << std::endl; LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sol_mdot)); LibmeshPetscCall(VecCopy(loc_vecs[0], sol_mdot)); LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &sol_p)); LibmeshPetscCall(VecCopy(loc_vecs[1], sol_p)); LibmeshPetscCall(VecDuplicate(_cmc_sys_Wij_rhs, &sol_Wij)); LibmeshPetscCall(VecCopy(loc_vecs[2], sol_Wij)); - if (_verbose_subchannel) - _console << "Getting individual field solutions from coupled solver." << std::endl; - - /// Assigning the solutions to arrays - PetscScalar * sol_p_array; - LibmeshPetscCall(VecGetArray(sol_p, &sol_p_array)); - PetscScalar * sol_Wij_array; - LibmeshPetscCall(VecGetArray(sol_Wij, &sol_Wij_array)); + V("Solution from coupled solver copied to solution vectors"); - /// Populating Mass flow + // mass flow LibmeshPetscCall(populateSolutionChan( sol_mdot, *_mdot_soln, first_node, last_node, _n_channels)); - /// Populating Pressure - for (unsigned int iz = last_node; iz > first_node - 1; iz--) + // pressure { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1); - PetscScalar value = sol_p_array[iz_ind * _n_channels + i_ch]; - _P_soln->set(node_in, value); - } - } - - /// Populating Crossflow - LibmeshPetscCall(populateDenseFromVector>( - sol_Wij, _Wij, first_node, last_node - 1, _n_gaps)); - - /// Populating Enthalpy - if (_monolithic_thermal_bool) - { - if (lag_block_thermal_solve) - { - KSP ksploc; - PC pc; - Vec sol; - LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol)); - LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc)); - LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat)); - LibmeshPetscCall(KSPGetPC(ksploc, &pc)); - LibmeshPetscCall(PCSetType(pc, PCJACOBI)); - LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit)); - LibmeshPetscCall(KSPSetFromOptions(ksploc)); - LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol)); - PetscScalar * xx; - LibmeshPetscCall(VecGetArray(sol, &xx)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) - { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); - auto h_out = xx[iz_ind * _n_channels + i_ch]; - if (h_out < 0) - { - mooseError(name(), - " : Calculation of negative Enthalpy h_out = : ", - h_out, - " Axial Level= : ", - iz); - } - _h_soln->set(node_out, h_out); - } - } - LibmeshPetscCall(KSPDestroy(&ksploc)); - LibmeshPetscCall(VecDestroy(&sol)); - } - else + PetscScalar * sol_p_array; + LibmeshPetscCall(VecGetArray(sol_p, &sol_p_array)); + for (unsigned int iz = last_node; iz > first_node - 1; iz--) { - Vec sol_h; - LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol_h)); - LibmeshPetscCall(VecCopy(loc_vecs[3], sol_h)); - PetscScalar * sol_h_array; - LibmeshPetscCall(VecGetArray(sol_h, &sol_h_array)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) + const auto iz_ind = iz - first_node; + for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); - auto h_out = sol_h_array[iz_ind * _n_channels + i_ch]; - if (h_out < 0) - { - mooseError(name(), - " : Calculation of negative Enthalpy h_out = : ", - h_out, - " Axial Level= : ", - iz); - } - _h_soln->set(node_out, h_out); - } + auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1); + PetscScalar value = sol_p_array[iz_ind * _n_channels + i_ch]; + _P_soln->set(node_in, value); } - LibmeshPetscCall(VecDestroy(&sol_h)); } + LibmeshPetscCall(VecRestoreArray(sol_p, &sol_p_array)); } - /// Populating sum_Wij + // crossflow dense + sumWij + correction factor + LibmeshPetscCall(populateDenseFromVector>( + sol_Wij, _Wij, first_node, last_node, _n_gaps)); + LibmeshPetscCall(MatMult(_mc_sumWij_mat, sol_Wij, _prod)); LibmeshPetscCall(populateSolutionChan( _prod, *_SumWij_soln, first_node, last_node, _n_channels)); - Vec sumWij_loc; - LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels)); - LibmeshPetscCall(populateVectorFromHandle( - _prod, *_SumWij_soln, first_node, last_node, _n_channels)); - LibmeshPetscCall(VecAbs(_prod)); LibmeshPetscCall(VecMax(_prod, NULL, &_max_sumWij_new)); - if (_verbose_subchannel) - _console << "Maximum estimated Wij new: " << _max_sumWij_new << std::endl; + V("Maximum estimated Wij new: " + std::to_string(_max_sumWij_new)); _correction_factor = _max_sumWij_new / _max_sumWij; - if (_verbose_subchannel) - _console << "Correction factor: " << _correction_factor << std::endl; - if (_verbose_subchannel) - _console << "Solutions assigned to MOOSE variables." << std::endl; + V("Correction factor: " + std::to_string(_correction_factor)); + V("Solutions assigned to MOOSE variables."); - /// Destroying arrays + // cleanup solution objects LibmeshPetscCall(VecDestroy(&x_nest)); LibmeshPetscCall(VecDestroy(&sol_mdot)); LibmeshPetscCall(VecDestroy(&sol_p)); LibmeshPetscCall(VecDestroy(&sol_Wij)); - LibmeshPetscCall(VecDestroy(&sumWij_loc)); - if (_verbose_subchannel) - _console << "Solutions destroyed." << std::endl; + V("Solutions destroyed."); PetscFunctionReturn(LIBMESH_PETSC_SUCCESS); } @@ -2648,8 +2513,14 @@ SubChannel1PhaseProblem::externalSolve() _dt = (isTransient() ? dt() : _one); _TR = isTransient(); initializeSolution(); - if (_verbose_subchannel) - _console << "Solution initialized" << std::endl; + // Small helper functions to reduce repetition + // Verbose print helper (no-op unless _verbose_subchannel is true) + auto V = [&](const std::string & s) + { + if (_verbose_subchannel) + _console << s << std::endl; + }; + V("Solution initialized"); Real P_error = 1.0; unsigned int P_it = 0; unsigned int P_it_max; @@ -2697,7 +2568,6 @@ SubChannel1PhaseProblem::externalSolve() if (_segregated_bool) { computeWijFromSolve(iblock); - if (_compute_power) { computeh(iblock); @@ -2708,38 +2578,21 @@ SubChannel1PhaseProblem::externalSolve() { LibmeshPetscCall(implicitPetscSolve(iblock)); computeWijPrime(iblock); - if (_verbose_subchannel) - _console << "Done with main solve." << std::endl; - if (_monolithic_thermal_bool) + V("Done with main solve."); + if (_compute_power) { - // Enthalpy is already solved from the monolithic solve + computeh(iblock); computeT(iblock); } - else - { - if (_verbose_subchannel) - _console << "Starting thermal solve." << std::endl; - if (_compute_power) - { - computeh(iblock); - computeT(iblock); - } - if (_verbose_subchannel) - _console << "Done with thermal solve." << std::endl; - } + V("Done with thermal solve."); } - if (_verbose_subchannel) - _console << "Start updating thermophysical properties." << std::endl; - + V("Start updating thermophysical properties."); if (_compute_density) computeRho(iblock); - if (_compute_viscosity) computeMu(iblock); - - if (_verbose_subchannel) - _console << "Done updating thermophysical properties." << std::endl; + V("Done updating thermophysical properties."); // We must do a global assembly to make sure data is parallel consistent before we do things // like compute L2 norms @@ -2759,16 +2612,21 @@ SubChannel1PhaseProblem::externalSolve() P_error = std::abs((P_L2norm_new_axial - P_L2norm_old_axial) / (P_L2norm_old_axial + _P_out + 1E-14)); _console << "P_error :" << P_error << std::endl; - if (_verbose_subchannel) - { - _console << "Iteration: " << P_it << std::endl; - _console << "Maximum iterations: " << P_it_max << std::endl; - } + V("Iteration: " + std::to_string(P_it)); + V("Maximum iterations: " + std::to_string(P_it_max)); } // update old crossflow matrix _Wij_old = _Wij; _console << "Finished executing subchannel solver\n"; + // set SumWij at the inlet equal to the one on the first axial level (for visualization purposes) + for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) + { + auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0); + auto * node_out = _subchannel_mesh.getChannelNode(i_ch, 1); + _SumWij_soln->set(node_in, (*_SumWij_soln)(node_out)); // kg/sec + } + /// Assigning temperature to the fuel pins if (_pin_mesh_exist) { diff --git a/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C b/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C index 9a37dcf5005e..9b69b9bc85bf 100644 --- a/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C +++ b/modules/subchannel/src/problems/TriSubChannel1PhaseProblem.C @@ -1476,44 +1476,8 @@ TriSubChannel1PhaseProblem::computeh(int iblock) LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_radial_heat_conduction_rhs)); LibmeshPetscCall(VecAXPY(_hc_sys_h_rhs, 1.0, _hc_sweep_enthalpy_rhs)); - if (_segregated_bool || (!_monolithic_thermal_bool)) - { - // Assembly the matrix system - KSP ksploc; - PC pc; - Vec sol; - LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol)); - LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc)); - LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat)); - LibmeshPetscCall(KSPGetPC(ksploc, &pc)); - LibmeshPetscCall(PCSetType(pc, PCJACOBI)); - LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit)); - LibmeshPetscCall(KSPSetOptionsPrefix(ksploc, "h_sys_")); - LibmeshPetscCall(KSPSetFromOptions(ksploc)); - LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol)); - // VecView(sol, PETSC_VIEWER_STDOUT_WORLD); - PetscScalar * xx; - LibmeshPetscCall(VecGetArray(sol, &xx)); - for (unsigned int iz = first_node; iz < last_node + 1; iz++) - { - auto iz_ind = iz - first_node; - for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++) - { - auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz); - auto h_out = xx[iz_ind * _n_channels + i_ch]; - if (h_out < 0) - { - mooseError(name(), - " : Calculation of negative Enthalpy h_out = : ", - h_out, - " Axial Level= : ", - iz); - } - _h_soln->set(node_out, h_out); - } - } - LibmeshPetscCall(KSPDestroy(&ksploc)); - LibmeshPetscCall(VecDestroy(&sol)); - } + // Use system to solve for and populate enthalpy + LibmeshPetscCall(this->solveAndPopulateEnthalpy( + _hc_sys_h_mat, _hc_sys_h_rhs, first_node, last_node, "h_sys_")); } } diff --git a/modules/subchannel/test/tests/positions/pin_positions.i b/modules/subchannel/test/tests/positions/pin_positions.i index 6ebec4312c28..70ac64c0b608 100644 --- a/modules/subchannel/test/tests/positions/pin_positions.i +++ b/modules/subchannel/test/tests/positions/pin_positions.i @@ -143,7 +143,6 @@ duct_inside = '${fparse duct_outside - 2 * duct_thickness}' implicit = false segregated = true staggered_pressure = false - monolithic_thermal = false # Tolerances P_tol = 1.0e-4 diff --git a/modules/subchannel/test/tests/problems/Lead-LBE-19pin/gold/test_LEAD-19pin_out.csv b/modules/subchannel/test/tests/problems/Lead-LBE-19pin/gold/test_LEAD-19pin_out.csv index 8be0abb01c78..718319f1c5b4 100644 --- a/modules/subchannel/test/tests/problems/Lead-LBE-19pin/gold/test_LEAD-19pin_out.csv +++ b/modules/subchannel/test/tests/problems/Lead-LBE-19pin/gold/test_LEAD-19pin_out.csv @@ -1,3 +1,3 @@ time,DP_SubchannelDelta,Mean_Temp,T1,T2,T3,T4,T5,T6,T7,T8,Total_power 0,0,0,0,0,0,0,0,0,0,0,0 -1,287402.54763993,725.29318657688,714.6413765836,714.46870030199,732.66221954555,735.30483779693,735.89369439068,735.89388897273,732.66276985976,714.64117590608,249999.99999999 +1,287404.44013129,725.29321735222,714.64116377983,714.46868862653,732.66305771052,735.30541398667,735.89356383531,735.8942811934,732.66258337294,714.64109329565,249999.99999999 diff --git a/modules/subchannel/test/tests/problems/Lead-LBE-19pin/test_LEAD-19pin.i b/modules/subchannel/test/tests/problems/Lead-LBE-19pin/test_LEAD-19pin.i index 8dba40c9d16c..90f47117d6d8 100644 --- a/modules/subchannel/test/tests/problems/Lead-LBE-19pin/test_LEAD-19pin.i +++ b/modules/subchannel/test/tests/problems/Lead-LBE-19pin/test_LEAD-19pin.i @@ -69,7 +69,6 @@ P_out = 1.0e5 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_multiapps = true verbose_subchannel = true interpolation_scheme = upwind diff --git a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_full_monolithic_out.csv b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_full_monolithic_out.csv deleted file mode 100644 index 45e723d5a653..000000000000 --- a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_full_monolithic_out.csv +++ /dev/null @@ -1,3 +0,0 @@ -time,DP_SubchannelDelta,Mean_Temp,T1,T2,T3,T4,T5,T6,T7,T8,Total_power,mdot-8 -0,0,0,0,0,0,0,0,0,0,0,0,0 -1,86121.575795134,660.37109543555,645.03559412795,663.41748118851,662.26769561659,661.67998595092,661.64388321406,661.64479519134,662.26818260387,645.03644609102,1000,0.1063894949395 diff --git a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_monolithic_out.csv b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_monolithic_out.csv index 7f92a4475218..bb93971b2e1a 100644 --- a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_monolithic_out.csv +++ b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/gold/test19_monolithic_out.csv @@ -1,3 +1,3 @@ time,DP_SubchannelDelta,Mean_Temp,T1,T2,T3,T4,T5,T6,T7,T8,Total_power,mdot-8 0,0,0,0,0,0,0,0,0,0,0,0,0 -1,86121.959155874,660.36796222015,645.0732145659,663.40571679344,662.26104026887,661.67742493403,661.64171240964,661.64160876006,662.26104529426,645.07320911146,1000,0.10639648579536 +1,86183.456095436,660.16180409665,660.12486187385,660.14641759136,660.18269103227,660.20959606091,660.24097838453,660.24106519599,660.18281326546,660.1249381626,1000,0.10754782191463 diff --git a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/test19_full_monolithic.i b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/test19_full_monolithic.i deleted file mode 100644 index b6e0d78ec7a9..000000000000 --- a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/test19_full_monolithic.i +++ /dev/null @@ -1,273 +0,0 @@ -T_in = 660 -mass_flux_in = '${fparse 1e+6 * 300.00 / 36000.*0.5}' -P_out = 2.0e5 # Pa - -[GlobalParams] - nrings = 3 - n_cells = 5 - flat_to_flat = 0.056 - heated_length = 0.5 - pitch = 0.012 -[] - -[TriSubChannelMesh] - [subchannel] - type = SCMTriSubChannelMeshGenerator - pin_diameter = 0.01 - dwire = 0.002 - hwire = 0.0833 - spacer_z = '0' - spacer_k = '5.0' - [] - [duct] - type = SCMTriDuctMeshGenerator - input = subchannel - [] -[] - -[AuxVariables] - [mdot] - block = subchannel - [] - [SumWij] - block = subchannel - [] - [P] - block = subchannel - [] - [DP] - block = subchannel - [] - [h] - block = subchannel - [] - [T] - block = subchannel - [] - [rho] - block = subchannel - [] - [S] - block = subchannel - [] - [w_perim] - block = subchannel - [] - [mu] - block = subchannel - [] - [q_prime] - block = subchannel - [] - [displacement] - block = subchannel - [] - [duct_heat_flux] - block = duct - [] - [Tduct] - block = duct - [] -[] - -[FluidProperties] - [sodium] - type = PBSodiumFluidProperties - [] -[] - -[Problem] - type = TriSubChannel1PhaseProblem - fp = sodium - n_blocks = 1 - P_out = 2.0e5 - CT = 1.0 - compute_density = true - compute_viscosity = true - compute_power = true - T_tol = 1.0e-6 - P_tol = 1.0e-6 - implicit = true - segregated = false - monolithic_thermal = true - - # Heat Transfer Correlations - duct_htc_correlation = 'gnielinski' -[] - -[ICs] - [S_IC] - type = SCMTriFlowAreaIC - variable = S - [] - - [w_perim_IC] - type = SCMTriWettedPerimIC - variable = w_perim - [] - - [q_prime_IC] - type = SCMTriPowerIC - variable = q_prime - power = 1000.0 # W - filename = "pin_power_profile19.txt" - [] - - [T_ic] - type = ConstantIC - variable = T - value = ${T_in} - [] - - [P_ic] - type = ConstantIC - variable = P - value = 0.0 - [] - - [DP_ic] - type = ConstantIC - variable = DP - value = 0.0 - [] - - [Viscosity_ic] - type = ViscosityIC - variable = mu - p = ${P_out} - T = T - fp = sodium - [] - - [rho_ic] - type = RhoFromPressureTemperatureIC - variable = rho - p = ${P_out} - T = T - fp = sodium - [] - - [h_ic] - type = SpecificEnthalpyFromPressureTemperatureIC - variable = h - p = ${P_out} - T = T - fp = sodium - [] - - [mdot_ic] - type = ConstantIC - variable = mdot - value = 0.0 - [] -[] - -[AuxKernels] - [T_in_bc] - type = ConstantAux - variable = T - boundary = inlet - value = ${T_in} - execute_on = 'timestep_begin' - [] - [mdot_in_bc] - type = SCMMassFlowRateAux - variable = mdot - boundary = inlet - area = S - mass_flux = ${mass_flux_in} - execute_on = 'timestep_begin' - [] -[] - -[Outputs] - exodus = true - csv = true -[] - -[Postprocessors] - [T1] - type = SubChannelPointValue - variable = T - index = 37 - execute_on = "timestep_end" - height = 0.5 - [] - [T2] - type = SubChannelPointValue - variable = T - index = 36 - execute_on = "timestep_end" - height = 0.5 - [] - [T3] - type = SubChannelPointValue - variable = T - index = 20 - execute_on = "timestep_end" - height = 0.5 - [] - [T4] - type = SubChannelPointValue - variable = T - index = 10 - execute_on = "timestep_end" - height = 0.5 - [] - [T5] - type = SubChannelPointValue - variable = T - index = 4 - execute_on = "timestep_end" - height = 0.5 - [] - [T6] - type = SubChannelPointValue - variable = T - index = 1 - execute_on = "timestep_end" - height = 0.5 - [] - [T7] - type = SubChannelPointValue - variable = T - index = 14 - execute_on = "timestep_end" - height = 0.5 - [] - [T8] - type = SubChannelPointValue - variable = T - index = 28 - execute_on = "timestep_end" - height = 0.5 - [] - ####### Assembly pressure drop - [DP_SubchannelDelta] - type = SubChannelDelta - variable = P - execute_on = 'TIMESTEP_END' - [] - ##### - [Mean_Temp] - type = SCMPlanarMean - variable = T - height = 2 - [] - [Total_power] - type = ElementIntegralVariablePostprocessor - variable = q_prime - block = subchannel - [] - [mdot-8] - type = SubChannelPointValue - variable = mdot - index = 28 - execute_on = 'TIMESTEP_END' - height = 0.5 - [] -[] - -[Executioner] - type = Steady -[] diff --git a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/tests b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/tests index 97a5715bee29..23e7f2db9e2d 100644 --- a/modules/subchannel/test/tests/problems/SFR/sodium-19pin/tests +++ b/modules/subchannel/test/tests/problems/SFR/sodium-19pin/tests @@ -29,14 +29,4 @@ recover = false requirement = 'The system will examine the subchannel monolithic solver for a 19pin liquid sodium cooled assembly' [] - [test_full_monolithic] - type = CSVDiff - input = test19_full_monolithic.i - csvdiff = test19_full_monolithic_out.csv - abs_zero = 1e-7 - rel_err = 5e-04 - valgrind = NONE - recover = false - requirement = 'The system will examine the subchannel full monolithic solver for a 19pin liquid sodium cooled assembly' - [] [] diff --git a/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_out.csv b/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_out.csv deleted file mode 100644 index 442a3e174147..000000000000 --- a/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_out.csv +++ /dev/null @@ -1,3 +0,0 @@ -time,T1,T2,T3,T4,T5,T6,report_mass_flux_inlet,total_pressure_drop -0,0,0,0,0,0,0,4722.2222222222,0 -1,387.60736447843,393.75458067795,387.42426636333,375.21721882344,368.58840419192,366.44585455997,4722.2222222222,29695.408166908 diff --git a/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_staggered_out.csv b/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_staggered_out.csv deleted file mode 100644 index 7125431a748e..000000000000 --- a/modules/subchannel/test/tests/problems/psbt/gold/psbt_full_monolithic_staggered_out.csv +++ /dev/null @@ -1,3 +0,0 @@ -time,T1,T2,T3,T4,T5,T6,report_mass_flux_inlet,total_pressure_drop -0,0,0,0,0,0,0,4722.2222222222,0 -1,387.64538423756,393.76103884862,387.38649800943,375.18939310852,368.59338680518,366.48192599777,4722.2222222222,29694.434544684 diff --git a/modules/subchannel/test/tests/problems/psbt/gold/psbt_monolithic_out.csv b/modules/subchannel/test/tests/problems/psbt/gold/psbt_monolithic_out.csv index f05bd3a1abc9..c11bff5d651d 100644 --- a/modules/subchannel/test/tests/problems/psbt/gold/psbt_monolithic_out.csv +++ b/modules/subchannel/test/tests/problems/psbt/gold/psbt_monolithic_out.csv @@ -1,3 +1,3 @@ time,T1,T2,T3,T4,T5,T6,report_mass_flux_inlet,total_pressure_drop 0,0,0,0,0,0,0,4722.2222222222,0 -1,387.60966732432,393.75694339233,387.42624035538,375.21852398553,368.59101187444,366.44702281277,4722.2222222222,29695.427851457 +1,387.50416729264,393.76233251255,387.51799410822,375.29284834476,368.59698496689,366.33664572405,4722.2222222222,29697.179314911 diff --git a/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic.i b/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic.i deleted file mode 100644 index e5bb7c3095ab..000000000000 --- a/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic.i +++ /dev/null @@ -1,208 +0,0 @@ -T_in = 359.15 -# [1e+6 kg/m^2-hour] turns into kg/m^2-sec -mass_flux_in = '${fparse 1e+6 * 17.00 / 3600.}' -P_out = 4.923e6 # Pa -pin_diameter = 0.00950 - -[QuadSubChannelMesh] - [sub_channel] - type = SCMQuadSubChannelMeshGenerator - nx = 6 - ny = 6 - n_cells = 10 - pitch = 0.0126 - pin_diameter = ${pin_diameter} - side_gap = 0.00095 - heated_length = 1.0 - spacer_z = '0.0' - spacer_k = '0.0' - [] - - [fuel_pins] - type = SCMQuadPinMeshGenerator - input = sub_channel - nx = 6 - ny = 6 - n_cells = 10 - pitch = 0.0126 - heated_length = 1.0 - [] -[] - -[FluidProperties] - [water] - type = Water97FluidProperties - [] -[] - -[SubChannel] - type = QuadSubChannel1PhaseProblem - fp = water - n_blocks = 1 - beta = 0.006 - CT = 2.6 - compute_density = true - compute_viscosity = true - compute_power = true - P_out = ${P_out} - verbose_subchannel = true - implicit = true - segregated = false - monolithic_thermal = true -[] - -[ICs] - [S_IC] - type = SCMQuadFlowAreaIC - variable = S - [] - - [w_perim_IC] - type = SCMQuadWettedPerimIC - variable = w_perim - [] - - [q_prime_IC] - type = SCMQuadPowerIC - variable = q_prime - power = 1.0e6 # W - filename = "power_profile.txt" #type in name of file that describes radial power profile - block = fuel_pins - [] - - [T_ic] - type = ConstantIC - variable = T - value = ${T_in} - [] - - [Dpin_ic] - type = ConstantIC - variable = Dpin - value = ${pin_diameter} - [] - - [P_ic] - type = ConstantIC - variable = P - value = 0.0 - [] - - [DP_ic] - type = ConstantIC - variable = DP - value = 0.0 - [] - - [Viscosity_ic] - type = ViscosityIC - variable = mu - p = ${P_out} - T = T - fp = water - [] - - [rho_ic] - type = RhoFromPressureTemperatureIC - variable = rho - p = ${P_out} - T = T - fp = water - [] - - [h_ic] - type = SpecificEnthalpyFromPressureTemperatureIC - variable = h - p = ${P_out} - T = T - fp = water - [] - - [mdot_ic] - type = ConstantIC - variable = mdot - value = 0.0 - [] -[] - -[AuxKernels] - [T_in_bc] - type = ConstantAux - variable = T - boundary = inlet - value = ${T_in} - execute_on = 'timestep_begin' - block = sub_channel - [] - [mdot_in_bc] - type = SCMMassFlowRateAux - variable = mdot - boundary = inlet - area = S - mass_flux = report_mass_flux_inlet - execute_on = 'timestep_begin' - block = sub_channel - [] -[] - -[Postprocessors] - [report_mass_flux_inlet] - type = Receiver - default = ${mass_flux_in} - [] - [total_pressure_drop] - type = SubChannelDelta - variable = P - execute_on = "timestep_end" - [] - [T1] - type = SubChannelPointValue - variable = T - index = 0 - execute_on = "timestep_end" - height = 1 - [] - [T2] - type = SubChannelPointValue - variable = T - index = 7 - execute_on = "timestep_end" - height = 1 - [] - [T3] - type = SubChannelPointValue - variable = T - index = 14 - execute_on = "timestep_end" - height = 1 - [] - [T4] - type = SubChannelPointValue - variable = T - index = 21 - execute_on = "timestep_end" - height = 1 - [] - [T5] - type = SubChannelPointValue - variable = T - index = 28 - execute_on = "timestep_end" - height = 1 - [] - [T6] - type = SubChannelPointValue - variable = T - index = 35 - execute_on = "timestep_end" - height = 1 - [] -[] - -[Outputs] - csv = true -[] - -[Executioner] - type = Steady -[] diff --git a/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic_staggered.i b/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic_staggered.i deleted file mode 100644 index d38fe59a3821..000000000000 --- a/modules/subchannel/test/tests/problems/psbt/psbt_full_monolithic_staggered.i +++ /dev/null @@ -1,209 +0,0 @@ -T_in = 359.15 -# [1e+6 kg/m^2-hour] turns into kg/m^2-sec -mass_flux_in = '${fparse 1e+6 * 17.00 / 3600.}' -P_out = 4.923e6 # Pa -pin_diameter = 0.00950 - -[QuadSubChannelMesh] - [sub_channel] - type = SCMQuadSubChannelMeshGenerator - nx = 6 - ny = 6 - n_cells = 10 - pitch = 0.0126 - pin_diameter = ${pin_diameter} - side_gap = 0.00095 - heated_length = 1.0 - spacer_z = '0.0' - spacer_k = '0.0' - [] - - [fuel_pins] - type = SCMQuadPinMeshGenerator - input = sub_channel - nx = 6 - ny = 6 - n_cells = 10 - pitch = 0.0126 - heated_length = 1.0 - [] -[] - -[FluidProperties] - [water] - type = Water97FluidProperties - [] -[] - -[SubChannel] - type = QuadSubChannel1PhaseProblem - fp = water - n_blocks = 1 - beta = 0.006 - CT = 2.6 - compute_density = true - compute_viscosity = true - compute_power = true - P_out = ${P_out} - verbose_subchannel = true - implicit = true - segregated = false - monolithic_thermal = true - staggered_pressure = true -[] - -[ICs] - [S_IC] - type = SCMQuadFlowAreaIC - variable = S - [] - - [w_perim_IC] - type = SCMQuadWettedPerimIC - variable = w_perim - [] - - [q_prime_IC] - type = SCMQuadPowerIC - variable = q_prime - power = 1.0e6 # W - filename = "power_profile.txt" #type in name of file that describes radial power profile - block = fuel_pins - [] - - [T_ic] - type = ConstantIC - variable = T - value = ${T_in} - [] - - [Dpin_ic] - type = ConstantIC - variable = Dpin - value = ${pin_diameter} - [] - - [P_ic] - type = ConstantIC - variable = P - value = 0.0 - [] - - [DP_ic] - type = ConstantIC - variable = DP - value = 0.0 - [] - - [Viscosity_ic] - type = ViscosityIC - variable = mu - p = ${P_out} - T = T - fp = water - [] - - [rho_ic] - type = RhoFromPressureTemperatureIC - variable = rho - p = ${P_out} - T = T - fp = water - [] - - [h_ic] - type = SpecificEnthalpyFromPressureTemperatureIC - variable = h - p = ${P_out} - T = T - fp = water - [] - - [mdot_ic] - type = ConstantIC - variable = mdot - value = 0.0 - [] -[] - -[AuxKernels] - [T_in_bc] - type = ConstantAux - variable = T - boundary = inlet - value = ${T_in} - execute_on = 'timestep_begin' - block = sub_channel - [] - [mdot_in_bc] - type = SCMMassFlowRateAux - variable = mdot - boundary = inlet - area = S - mass_flux = report_mass_flux_inlet - execute_on = 'timestep_begin' - block = sub_channel - [] -[] - -[Postprocessors] - [report_mass_flux_inlet] - type = Receiver - default = ${mass_flux_in} - [] - [total_pressure_drop] - type = SubChannelDelta - variable = P - execute_on = "timestep_end" - [] - [T1] - type = SubChannelPointValue - variable = T - index = 0 - execute_on = "timestep_end" - height = 1 - [] - [T2] - type = SubChannelPointValue - variable = T - index = 7 - execute_on = "timestep_end" - height = 1 - [] - [T3] - type = SubChannelPointValue - variable = T - index = 14 - execute_on = "timestep_end" - height = 1 - [] - [T4] - type = SubChannelPointValue - variable = T - index = 21 - execute_on = "timestep_end" - height = 1 - [] - [T5] - type = SubChannelPointValue - variable = T - index = 28 - execute_on = "timestep_end" - height = 1 - [] - [T6] - type = SubChannelPointValue - variable = T - index = 35 - execute_on = "timestep_end" - height = 1 - [] -[] - -[Outputs] - csv = true -[] - -[Executioner] - type = Steady -[] diff --git a/modules/subchannel/test/tests/problems/psbt/tests b/modules/subchannel/test/tests/problems/psbt/tests index 4110de5f2a66..7d288772161a 100644 --- a/modules/subchannel/test/tests/problems/psbt/tests +++ b/modules/subchannel/test/tests/problems/psbt/tests @@ -56,28 +56,4 @@ max_threads = 1 requirement = 'The system shall be able to solve a PSBT-type case with the subchannel solver using a monolithic algorithm.' [] - [psbt_regression_test_full_monolithic] - type = CSVDiff - input = psbt_full_monolithic.i - csvdiff = psbt_full_monolithic_out.csv - capabilities = 'method!=dbg' - valgrind = NONE - recover = false - abs_zero = 1e-6 - rel_err = 1e-4 - max_threads = 1 - requirement = 'The system shall be able to solve a PSBT-type case with the subchannel solver using a full-monolithic algorithm.' - [] - [psbt_regression_test_full_monolithic_staggered] - type = CSVDiff - input = psbt_full_monolithic_staggered.i - csvdiff = psbt_full_monolithic_staggered_out.csv - capabilities = 'method!=dbg' - valgrind = NONE - recover = false - abs_zero = 1e-6 - rel_err = 6e-5 - max_threads = 1 - requirement = 'The system shall be able to solve a PSBT-type case with the subchannel solver using a full-monolithic algorithm and a staggered pressure formulation.' - [] [] diff --git a/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage70.i b/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage70.i index e88ca4bffc34..e35728825b62 100644 --- a/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage70.i +++ b/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage70.i @@ -77,7 +77,6 @@ P_out = 101325 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false interpolation_scheme = 'upwind' [] diff --git a/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage90.i b/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage90.i index e6012365aa20..33a1addc1570 100644 --- a/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage90.i +++ b/modules/subchannel/validation/Blockage/PNNL_7x7/7X7blockage90.i @@ -77,7 +77,6 @@ P_out = 101325 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false interpolation_scheme = central_difference [] diff --git a/modules/subchannel/validation/Blockage/THORS/gold/FFM-3A_out.csv b/modules/subchannel/validation/Blockage/THORS/gold/FFM-3A_out.csv index 815a8a72d87b..e9fd75a2ee01 100644 --- a/modules/subchannel/validation/Blockage/THORS/gold/FFM-3A_out.csv +++ b/modules/subchannel/validation/Blockage/THORS/gold/FFM-3A_out.csv @@ -1,3 +1,3 @@ time,1,2,3,4,5,6,7,8 0,0,0,0,0,0,0,0,0 -1,790.30187042744,791.90405425724,813.38861804593,827.07682228701,841.84520534379,841.84422101146,813.38764658885,790.30212137055 +1,790.24655874655,791.81333809364,813.25122774133,827.12133000455,842.50021608762,842.50087363271,813.25146050529,790.2469707299 diff --git a/modules/subchannel/validation/ORNL_19_pin/ORNL_19.i b/modules/subchannel/validation/ORNL_19_pin/ORNL_19.i index 31096ca585d8..fb8455144cf0 100644 --- a/modules/subchannel/validation/ORNL_19_pin/ORNL_19.i +++ b/modules/subchannel/validation/ORNL_19_pin/ORNL_19.i @@ -66,7 +66,6 @@ P_out = 2.0e5 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = true verbose_multiapps = true verbose_subchannel = true interpolation_scheme = upwind diff --git a/modules/subchannel/validation/ORNL_19_pin/test_ORNL_19.i b/modules/subchannel/validation/ORNL_19_pin/test_ORNL_19.i index 7b68330c1f85..d7b26c61b79b 100644 --- a/modules/subchannel/validation/ORNL_19_pin/test_ORNL_19.i +++ b/modules/subchannel/validation/ORNL_19_pin/test_ORNL_19.i @@ -74,7 +74,6 @@ P_out = 2.0e5 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_multiapps = true verbose_subchannel = false [] diff --git a/modules/subchannel/validation/PNNL_12_pin/steady_state/2X6_ss.i b/modules/subchannel/validation/PNNL_12_pin/steady_state/2X6_ss.i index 75587a8a650d..8e5ed4bd3a8b 100644 --- a/modules/subchannel/validation/PNNL_12_pin/steady_state/2X6_ss.i +++ b/modules/subchannel/validation/PNNL_12_pin/steady_state/2X6_ss.i @@ -88,7 +88,6 @@ P_out = 101325 # Pa P_out = ${P_out} implicit = true segregated = false - monolithic_thermal = false [] [ICs] diff --git a/modules/subchannel/validation/PNNL_12_pin/steady_state/gold/2X6_ss_out.csv b/modules/subchannel/validation/PNNL_12_pin/steady_state/gold/2X6_ss_out.csv index f9f6ad681283..f35ace3923a4 100644 --- a/modules/subchannel/validation/PNNL_12_pin/steady_state/gold/2X6_ss_out.csv +++ b/modules/subchannel/validation/PNNL_12_pin/steady_state/gold/2X6_ss_out.csv @@ -1,3 +1,3 @@ time,mdot10,mdot11,mdot12,mdot13,mdot7,mdot8,mdot9 0,0,0,0,0,0,0,0 -1,0.014962533415434,0.016611852105515,0.016423871332285,0.0095423970991296,0.0088959670195205,0.014246404016697,0.013641513521631 +1,0.014962354677177,0.016611382231055,0.016423580385621,0.0095424939657114,0.0088959970433663,0.01424633349384,0.013641594900661 diff --git a/modules/subchannel/validation/Toshiba_37_pin/gold/toshiba_37_pin_out.csv b/modules/subchannel/validation/Toshiba_37_pin/gold/toshiba_37_pin_out.csv index 8b104f4a1be0..5a7994582cfc 100644 --- a/modules/subchannel/validation/Toshiba_37_pin/gold/toshiba_37_pin_out.csv +++ b/modules/subchannel/validation/Toshiba_37_pin/gold/toshiba_37_pin_out.csv @@ -1,3 +1,3 @@ time,DP_SubchannelDelta,T_Planar_Mean 0,0,0 -1,9286.5384435489,707.4469412441 +1,9286.6619132485,707.44050544809 diff --git a/modules/subchannel/validation/Toshiba_37_pin/toshiba_37_pin.i b/modules/subchannel/validation/Toshiba_37_pin/toshiba_37_pin.i index 8291e73dbfe3..281f3fa14813 100644 --- a/modules/subchannel/validation/Toshiba_37_pin/toshiba_37_pin.i +++ b/modules/subchannel/validation/Toshiba_37_pin/toshiba_37_pin.i @@ -65,7 +65,6 @@ P_out = 2.0e5 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_multiapps = true verbose_subchannel = false [] diff --git a/modules/subchannel/validation/psbt/psbt_ss/gold/psbt_out.csv b/modules/subchannel/validation/psbt/psbt_ss/gold/psbt_out.csv index dd456185febc..9b2e74c7ca7e 100644 --- a/modules/subchannel/validation/psbt/psbt_ss/gold/psbt_out.csv +++ b/modules/subchannel/validation/psbt/psbt_ss/gold/psbt_out.csv @@ -1,3 +1,3 @@ time,PinTemp,T1,T2,T3,T4,T5,T6,total_pressure_drop 0,0,0,0,0,0,0,0,0 -1,612.02990188326,567.99855284885,568.5500625951,564.70705442123,558.19402276689,552.63785457705,549.84218031668,255377.02437708 +1,612.02925913972,568.66454103856,568.50025902784,564.3993713214,557.86071337793,552.59189010377,550.56679424957,255376.14242658 diff --git a/modules/subchannel/validation/psbt/psbt_ss/psbt.i b/modules/subchannel/validation/psbt/psbt_ss/psbt.i index 30d56d469b76..73f0d3d560a2 100644 --- a/modules/subchannel/validation/psbt/psbt_ss/psbt.i +++ b/modules/subchannel/validation/psbt/psbt_ss/psbt.i @@ -89,7 +89,6 @@ P_out = 14.72e6 # Pa implicit = true segregated = false staggered_pressure = false - monolithic_thermal = false verbose_subchannel = true interpolation_scheme = exponential deformation = true # this flag allows the re-calculation of subchannel geometric parameters based on the dpin value