https://mooseframework.inl.gov
SubChannel1PhaseProblem.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #pragma once
11 #include "ExternalProblem.h"
12 #include "PostprocessorInterface.h"
13 #include "SubChannelApp.h"
14 #include "QuadSubChannelMesh.h"
15 #include "SolutionHandle.h"
17 #include <petscdm.h>
18 #include <petscdmda.h>
19 #include <petscksp.h>
20 #include <petscsys.h>
21 #include <petscvec.h>
22 #include <petscsnes.h>
23 
28 {
29 public:
31  virtual ~SubChannel1PhaseProblem();
32 
33  virtual void externalSolve() override;
34  virtual void syncSolutions(Direction direction) override;
35  virtual bool solverSystemConverged(const unsigned int) override;
36  virtual void initialSetup() override;
37 
39  virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) = 0;
40 
41 protected:
43  {
44  int i_ch;
47 
49  virtual Real computeFrictionFactor(FrictionStruct friction_args) = 0;
51  void computeWijFromSolve(int iblock);
53  void computeSumWij(int iblock);
55  void computeMdot(int iblock);
57  void computeWijPrime(int iblock);
59  virtual Real computeBeta(unsigned int i_gap, unsigned int iz) = 0;
61  void computeDP(int iblock);
63  void computeP(int iblock);
65  virtual void computeh(int iblock) = 0;
67  void computeT(int iblock);
69  void computeRho(int iblock);
71  void computeMu(int iblock);
73  void computeWijResidual(int iblock);
75  Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz);
79  PetscErrorCode petscSnesSolver(int iblock,
80  const libMesh::DenseVector<Real> & solution,
83  friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void * ctx);
84 
86  PetscErrorCode implicitPetscSolve(int iblock);
87 
89  virtual void initializeSolution() = 0;
90 
92  PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
93  PetscScalar
94  computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
95 
96  PetscErrorCode cleanUp();
99  unsigned int _n_blocks;
104  const Real _g_grav;
105  const Real & _kij;
106  unsigned int _n_cells;
107  unsigned int _n_gaps;
108  unsigned int _n_pins;
109  unsigned int _n_channels;
110  unsigned int _block_size;
113  std::vector<Real> _z_grid;
118  const bool _compute_density;
120  const bool _compute_viscosity;
122  const bool _compute_power;
124  const bool _pin_mesh_exist;
126  const bool _duct_mesh_exist;
134  const Real & _CT;
136  const Real & _P_tol;
138  const Real & _T_tol;
140  const int & _T_maxit;
142  const PetscReal & _rtol;
144  const PetscReal & _atol;
146  const PetscReal & _dtol;
148  const PetscInt & _maxit;
152  const bool _implicit_bool;
156  const bool _segregated_bool;
162  const bool _deformation;
163 
166  std::unique_ptr<SolutionHandle> _mdot_soln;
167  std::unique_ptr<SolutionHandle> _SumWij_soln;
168  std::unique_ptr<SolutionHandle> _P_soln;
169  std::unique_ptr<SolutionHandle> _DP_soln;
170  std::unique_ptr<SolutionHandle> _h_soln;
171  std::unique_ptr<SolutionHandle> _T_soln;
172  std::unique_ptr<SolutionHandle> _Tpin_soln;
173  std::unique_ptr<SolutionHandle> _Dpin_soln;
174  std::unique_ptr<SolutionHandle> _rho_soln;
175  std::unique_ptr<SolutionHandle> _mu_soln;
176  std::unique_ptr<SolutionHandle> _S_flow_soln;
177  std::unique_ptr<SolutionHandle> _w_perim_soln;
178  std::unique_ptr<SolutionHandle> _q_prime_soln;
179  std::unique_ptr<SolutionHandle> _q_prime_duct_soln; // Only used for ducted assemblies
180  std::unique_ptr<SolutionHandle> _Tduct_soln; // Only used for ducted assemblies
181  std::unique_ptr<SolutionHandle> _displacement_soln;
182 
184  inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
185  {
187  LibmeshPetscCall(VecCreate(PETSC_COMM_WORLD, &v));
188  LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
189  LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
190  LibmeshPetscCall(VecSetFromOptions(v));
191  LibmeshPetscCall(VecZeroEntries(v));
192  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
193  }
194 
195  inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
196  {
198  LibmeshPetscCall(MatCreate(PETSC_COMM_WORLD, &M));
199  LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
200  LibmeshPetscCall(MatSetFromOptions(M));
201  LibmeshPetscCall(MatSetUp(M));
202  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
203  }
204 
205  template <class T>
206  PetscErrorCode populateVectorFromDense(Vec & x,
207  const T & solution,
208  const unsigned int first_axial_level,
209  const unsigned int last_axial_level,
210  const unsigned int cross_dimension);
211  template <class T>
212  PetscErrorCode populateDenseFromVector(const Vec & x,
213  T & solution,
214  const unsigned int first_axial_level,
215  const unsigned int last_axial_level,
216  const unsigned int cross_dimension);
217  template <class T>
218  PetscErrorCode populateVectorFromHandle(Vec & x,
219  const T & solution,
220  const unsigned int first_axial_level,
221  const unsigned int last_axial_level,
222  const unsigned int cross_dimension);
223  template <class T>
224  PetscErrorCode populateSolutionChan(const Vec & x,
225  T & solution,
226  const unsigned int first_axial_level,
227  const unsigned int last_axial_level,
228  const unsigned int cross_dimension);
229 
230  template <class T>
231  PetscErrorCode populateSolutionGap(const Vec & x,
232  T & solution,
233  const unsigned int first_axial_level,
234  const unsigned int last_axial_level,
235  const unsigned int cross_dimension);
236 
241  Vec _Wij_vec;
242  Vec _prod;
243  Vec _prodp;
249 
275 
292 
308 
310  PetscScalar _added_K = 0.0;
311  PetscScalar _added_K_old = 1000.0;
312  PetscScalar _max_sumWij;
313  PetscScalar _max_sumWij_new;
314  PetscScalar _correction_factor = 1.0;
315 
316 public:
317  static InputParameters validParams();
318 };
319 
320 template <class T>
321 PetscErrorCode
323  T & loc_solution,
324  const unsigned int first_axial_level,
325  const unsigned int last_axial_level,
326  const unsigned int cross_dimension)
327 {
328  PetscScalar * xx;
329 
331  LibmeshPetscCall(VecGetArray(x, &xx));
332  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
333  {
334  unsigned int iz_ind = iz - first_axial_level;
335  for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
336  {
337  loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
338  }
339  }
340  LibmeshPetscCall(VecRestoreArray(x, &xx));
341 
342  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
343 }
344 
345 template <class T>
346 PetscErrorCode
348  const T & loc_solution,
349  const unsigned int first_axial_level,
350  const unsigned int last_axial_level,
351  const unsigned int cross_dimension)
352 {
353  PetscScalar * xx;
354 
356  LibmeshPetscCall(VecGetArray(x, &xx));
357  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
358  {
359  unsigned int iz_ind = iz - first_axial_level;
360  for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
361  {
362  auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
363  xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
364  }
365  }
366  LibmeshPetscCall(VecRestoreArray(x, &xx));
367 
368  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
369 }
370 
371 template <class T>
372 PetscErrorCode
374  const T & loc_solution,
375  const unsigned int first_axial_level,
376  const unsigned int last_axial_level,
377  const unsigned int cross_dimension)
378 {
379  PetscScalar * xx;
381  LibmeshPetscCall(VecGetArray(x, &xx));
382  for (unsigned int iz = first_axial_level; iz < last_axial_level; iz++)
383  {
384  unsigned int iz_ind = iz - first_axial_level;
385  for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
386  {
387  xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
388  }
389  }
390  LibmeshPetscCall(VecRestoreArray(x, &xx));
391  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
392 }
393 
394 template <class T>
395 PetscErrorCode
397  T & loc_solution,
398  const unsigned int first_axial_level,
399  const unsigned int last_axial_level,
400  const unsigned int cross_dimension)
401 {
402  PetscScalar * xx;
404  LibmeshPetscCall(VecGetArray(x, &xx));
405  Node * loc_node;
406  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
407  {
408  unsigned int iz_ind = iz - first_axial_level;
409  for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
410  {
411  loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
412  loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
413  }
414  }
415  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
416 }
417 
418 template <class T>
419 PetscErrorCode
421  T & loc_solution,
422  const unsigned int first_axial_level,
423  const unsigned int last_axial_level,
424  const unsigned int cross_dimension)
425 {
426  PetscScalar * xx;
428  LibmeshPetscCall(VecGetArray(x, &xx));
429  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
430  {
431  unsigned int iz_ind = iz - first_axial_level;
432  for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
433  {
434  loc_solution(iz * cross_dimension + i_l) = xx[iz_ind * cross_dimension + i_l];
435  }
436  }
437  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
438 }
Mat _amc_sys_mdot_mat
Axial momentum system matrix.
const bool _pin_mesh_exist
Flag that informs if there is a pin mesh or not.
void computeRho(int iblock)
Computes Density per channel for block iblock.
Vec _amc_gravity_rhs
Axial momentum conservation - buoyancy force No implicit matrix.
std::unique_ptr< SolutionHandle > _Tduct_soln
void computeSumWij(int iblock)
Computes net diversion crossflow per channel for block iblock.
unsigned int _n_blocks
number of axial blocks
const bool _compute_power
Flag that informs if we need to solve the Enthalpy/Temperature equations or not.
std::unique_ptr< SolutionHandle > _T_soln
std::unique_ptr< SolutionHandle > _h_soln
const PetscReal & _dtol
The divergence tolerance for the ksp linear solver.
virtual Real computeBeta(unsigned int i_gap, unsigned int iz)=0
Computes turbulent mixing coefficient.
const bool _monolithic_thermal_bool
Thermal monolithic bool.
std::unique_ptr< SolutionHandle > _q_prime_duct_soln
Mat _cmc_friction_force_mat
Cross momentum conservation - friction force.
void computeP(int iblock)
Computes Pressure per channel for block iblock.
Real _TR
Flag that activates or deactivates the transient parts of the equations we solve by multiplication...
libMesh::DenseMatrix< Real > _Wij_old
const PostprocessorValue & _P_out
Outlet Pressure.
void computeT(int iblock)
Computes Temperature per channel for block iblock.
void computeMdot(int iblock)
Computes mass flow per channel for block iblock.
std::unique_ptr< SolutionHandle > _S_flow_soln
const bool _compute_density
Flag that activates or deactivates the calculation of density.
PetscErrorCode populateSolutionChan(const Vec &x, T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
Mat _amc_friction_force_mat
Axial momentum conservation - friction force.
PetscScalar _added_K
Added resistances for monolithic convergence.
PetscErrorCode createPetscVector(Vec &v, PetscInt n)
Petsc Functions.
PetscFunctionBegin
const bool _duct_mesh_exist
Flag that informs if there is a pin mesh or not.
static InputParameters validParams()
virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz)=0
Computes added heat for channel i_ch and cell iz.
const Real & _T_tol
Convergence tolerance for the temperature loop in external solve.
void computeMu(int iblock)
Computes Viscosity per channel for block iblock.
PetscScalar computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet=0.0)
const bool _segregated_bool
Segregated solve.
Mat _cmc_sys_Wij_mat
Lateral momentum system matrix.
bool _converged
Variable that informs whether we exited external solve with a converged solution or not...
const bool _staggered_pressure_bool
Flag to define the usage of staggered or collocated pressure.
virtual void initializeSolution()=0
Function to initialize the solution & geometry fields.
Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz)
Function that computes the heat flux added by the duct.
Mat _mc_sumWij_mat
Matrices and vectors to be used in implicit assembly Mass conservation Mass conservation - sum of cro...
PetscErrorCode populateVectorFromDense(Vec &x, const T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
const Real & _CT
Turbulent modeling parameter used in axial momentum equation.
virtual Node * getChannelNode(unsigned int i_chan, unsigned iz) const =0
Get the subchannel mesh node for a given channel index and elevation index.
virtual void computeh(int iblock)=0
Computes Enthalpy per channel for block iblock.
std::unique_ptr< SolutionHandle > _rho_soln
Mat _amc_turbulent_cross_flows_mat
Mass conservation - density time derivative No implicit matrix.
PetscErrorCode populateVectorFromHandle(Vec &x, const T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
Mat _amc_advective_derivative_mat
Axial momentum conservation - advective (Eulerian) derivative.
Vec _hc_added_heat_rhs
Enthalpy conservation - source and sink.
std::vector< Real > _z_grid
axial location of nodes
const int & _T_maxit
Maximum iterations for the inner temperature loop.
Mat _amc_pressure_force_mat
Axial momentum conservation - pressure force.
void computeWijResidual(int iblock)
Computes Residual Matrix based on the lateral momentum conservation equation for block iblock...
const std::vector< double > x
Real f(Real x)
Test function for Brents method.
std::unique_ptr< SolutionHandle > _q_prime_soln
virtual void syncSolutions(Direction direction) override
PetscErrorCode populateSolutionGap(const Vec &x, T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
const PetscReal & _atol
The absolute convergence tolerance for the ksp linear solver.
PetscScalar computeInterpolationCoefficients(PetscScalar Peclet=0.0)
Functions that computes the interpolation scheme given the Peclet number.
Real PostprocessorValue
Real root(std::function< Real(Real)> const &f, Real x1, Real x2, Real tol=1.0e-12)
Finds the root of a function using Brent&#39;s method.
Definition: BrentsMethod.C:66
Mat _hc_cross_derivative_mat
Enthalpy conservation - cross flux derivative.
Mat _mc_axial_convection_mat
Mass conservation - axial convection.
Common class for single phase fluid properties.
auto Peclet(const T1 &volume_fraction, const T2 &cp, const T3 &rho, const T4 &vel, const T5 &D_h, const T6 &k)
Compute Peclet number.
Definition: Numerics.h:153
virtual Real computeFrictionFactor(FrictionStruct friction_args)=0
Returns friction factor.
const PetscInt & _maxit
The maximum number of iterations to use for the ksp linear solver.
PetscErrorCode createPetscMatrix(Mat &M, PetscInt n, PetscInt m)
friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void *ctx)
This is the residual Vector function in a form compatible with the SNES PETC solvers.
const MooseEnum _interpolation_scheme
The interpolation method used in constructing the systems.
std::unique_ptr< SolutionHandle > _Dpin_soln
Base class for the 1-phase steady-state/transient subchannel solver.
libMesh::DenseMatrix< Real > _Wij_residual_matrix
std::unique_ptr< SolutionHandle > _mdot_soln
PetscErrorCode populateDenseFromVector(const Vec &x, T &solution, const unsigned int first_axial_level, const unsigned int last_axial_level, const unsigned int cross_dimension)
const PetscReal & _rtol
The relative convergence tolerance, (relative decrease) for the ksp linear solver.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unique_ptr< SolutionHandle > _displacement_soln
static const std::string v
Definition: NS.h:84
std::unique_ptr< SolutionHandle > _DP_soln
const bool _compute_viscosity
Flag that activates or deactivates the calculation of viscosity.
libMesh::DenseMatrix< Real > & _Wij
virtual void externalSolve() override
void computeWijPrime(int iblock)
Computes turbulent crossflow per gap for block iblock.
std::unique_ptr< SolutionHandle > _SumWij_soln
PetscErrorCode petscSnesSolver(int iblock, const libMesh::DenseVector< Real > &solution, libMesh::DenseVector< Real > &root)
Computes solution of nonlinear equation using snes and provided a residual in a formFunction.
libMesh::DenseVector< Real > residualFunction(int iblock, libMesh::DenseVector< Real > solution)
Computes Residual Vector based on the lateral momentum conservation equation for block iblock & updat...
const bool _verbose_subchannel
Boolean to printout information related to subchannel solve.
Mat _amc_cross_derivative_mat
Axial momentum conservation - cross flux derivative.
const Real & _P_tol
Convergence tolerance for the pressure loop in external solve.
virtual bool solverSystemConverged(const unsigned int) override
struct SubChannel1PhaseProblem::FrictionStruct _friction_args
Mat _hc_advective_derivative_mat
Enthalpy conservation - advective (Eulerian) derivative;.
Base class for subchannel meshes.
void computeWijFromSolve(int iblock)
Computes diversion crossflow per gap for block iblock.
const bool _implicit_bool
Flag to define the usage of a implicit or explicit solution.
Mat _amc_time_derivative_mat
Axial momentum conservation - time derivative.
Mat _hc_time_derivative_mat
Enthalpy Enthalpy conservation - time derivative.
PetscErrorCode implicitPetscSolve(int iblock)
Computes implicit solve using PetSc.
const bool _deformation
Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement...
void computeDP(int iblock)
Computes Pressure Drop per channel for block iblock.
SubChannel1PhaseProblem(const InputParameters &params)
Mat _cmc_time_derivative_mat
Cross momentum Cross momentum conservation - time derivative.
libMesh::DenseMatrix< Real > _WijPrime
virtual void initialSetup() override
PetscFunctionReturn(LIBMESH_PETSC_SUCCESS)
Mat _cmc_pressure_force_mat
Cross momentum conservation - pressure force.
std::unique_ptr< SolutionHandle > _P_soln
std::unique_ptr< SolutionHandle > _w_perim_soln
Mat _cmc_advective_derivative_mat
Cross momentum conservation - advective (Eulerian) derivative.
const SinglePhaseFluidProperties * _fp
Solutions handles and link to TH tables properties.
Mat _hc_sys_h_mat
System matrices.
std::unique_ptr< SolutionHandle > _mu_soln
std::unique_ptr< SolutionHandle > _Tpin_soln