diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 3544f028267..e21b2a559cc 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -135,7 +135,8 @@ class CConfig { Hold_GridFixed, /*!< \brief Flag hold fixed some part of the mesh during the deformation. */ Axisymmetric, /*!< \brief Flag for axisymmetric calculations */ Enable_Cuda, /*!< \brief Flag for switching GPU computing*/ - Integrated_HeatFlux; /*!< \brief Flag for heat flux BC whether it deals with integrated values.*/ + Integrated_HeatFlux, /*!< \brief Flag for heat flux BC whether it deals with integrated values.*/ + Pressure_Based; /*!< \brief Flag to check if we are using a pressure-based system.*/ su2double Buffet_k; /*!< \brief Sharpness coefficient for buffet sensor.*/ su2double Buffet_lambda; /*!< \brief Offset parameter for buffet sensor.*/ su2double Damp_Engine_Inflow; /*!< \brief Damping factor for the engine inlet. */ @@ -525,6 +526,8 @@ class CConfig { Kind_Gradient_Method_Recon, /*!< \brief Numerical method for computation of spatial gradients used for upwind reconstruction. */ Kind_Deform_Linear_Solver, /*!< Numerical method to deform the grid */ Kind_Deform_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ + Kind_Poisson_Linear_Solver, /*!< \brief Numerical solver for the poisson equation. */ + Kind_Poisson_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver of the poisson equation. */ Kind_Linear_Solver, /*!< \brief Numerical solver for the implicit scheme. */ Kind_Linear_Solver_Prec, /*!< \brief Preconditioner of the linear solver. */ Kind_DiscAdj_Linear_Solver, /*!< \brief Linear solver for the discrete adjoint system. */ @@ -587,6 +590,9 @@ class CConfig { Kind_Upwind_Heat, /*!< \brief Upwind scheme for the heat transfer model. */ Kind_Upwind_Template; /*!< \brief Upwind scheme for the template model. */ + PBITER Kind_PBIter; /*< \brief Kind of pressure-based algorithm that is used. */ + INCOMP_SYSTEM Kind_Incomp_System; /*< \brief Kind of incompressible solver. */ + bool MUSCL, /*!< \brief MUSCL scheme (for the runtime eq. system). */ MUSCL_Flow, /*!< \brief MUSCL scheme for the flow equations.*/ MUSCL_Turb, /*!< \brief MUSCL scheme for the turbulence equations.*/ @@ -638,8 +644,10 @@ class CConfig { bool InletUseNormal; /*!< \brief Flag for whether to use the local normal as the flow direction for a pressure inlet. */ su2double Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ + su2double Poisson_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the poisson equation. */ su2double Linear_Solver_Smoother_Relaxation; /*!< \brief Relaxation factor for iterative linear smoothers. */ unsigned long Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ + unsigned long Poisson_Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the poisson solver*/ unsigned long Deform_Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ unsigned long Linear_Solver_Restart_Frequency; /*!< \brief Restart frequency of the linear solver for the implicit formulation. */ unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */ @@ -655,6 +663,14 @@ class CConfig { su2double SemiSpan; /*!< \brief Wing Semi span. */ su2double MSW_Alpha; /*!< \brief Coefficient for blending states in the MSW scheme. */ su2double Roe_Kappa; /*!< \brief Relaxation of the Roe scheme. */ + + struct CSIMPLE_Options { + su2double Transient_Term_Removal_Factor; /*!< \brief Coefficient for removing the transient term from the momentum coefficient. */ + su2double Relaxation_Factor_Pressure; /*!< \brief Relaxation coefficient of the pressure corrections in the SIMPLE solver. */ + bool AutomaticRelaxationFactors; /*!< \brief option for automatically computing relaxation factors for flow corrections in SIMPLE. */ + unsigned short nCorrections_PISO; /*!< \brief Number of corrections used in PISO algorithm. */ + } SIMPLE_Options; + su2double Relaxation_Factor_Adjoint; /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */ su2double Relaxation_Factor_CHT; /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */ su2double EntropyFix_Coeff; /*!< \brief Entropy fix coefficient. */ @@ -4006,6 +4022,18 @@ class CConfig { */ ENUM_REGIME GetKind_Regime(void) const { return Kind_Regime; } + /*! + * \brief Kind of incompressible solver formulation. + * \return Kind of incompressible solver. + */ + INCOMP_SYSTEM GetKind_Incomp_System(void) const { return Kind_Incomp_System; } + + /*! + * \brief Kind of iteration used for pressure based iterations. + * \return Kind of iteration used for pressure based iterations. + */ + PBITER GetKind_PBIter(void) const { return Kind_PBIter; } + /*! * \brief Governing equations of the flow (it can be different from the run time equation). * \param[in] val_zone - Zone where the soler is applied. @@ -4356,12 +4384,24 @@ class CConfig { */ unsigned short GetKind_Linear_Solver_Prec(void) const { return Kind_Linear_Solver_Prec; } + /*! + * \brief Get the kind of preconditioner for the linear solver of the poisson problem. + * \return Numerical preconditioner for poisson equation (solving the linear system). + */ + unsigned short GetKind_Poisson_Linear_Solver_Prec(void) const { return Kind_Poisson_Linear_Solver_Prec; } + /*! * \brief Get the kind of solver for the implicit solver. * \return Numerical solver for implicit formulation (solving the linear system). */ unsigned short GetKind_Deform_Linear_Solver(void) const { return Kind_Deform_Linear_Solver; } + /*! + * \brief Get the kind of solver for the poisson equation. + * \return Numerical solver for poisson equation (solving the linear system). + */ + unsigned short GetKind_Poisson_Linear_Solver(void) const { return Kind_Poisson_Linear_Solver; } + /*! * \brief Get min error of the linear solver for the implicit formulation. * \return Min error of the linear solver for the implicit formulation. @@ -4374,12 +4414,24 @@ class CConfig { */ su2double GetDeform_Linear_Solver_Error(void) const { return Deform_Linear_Solver_Error; } + /*! + * \brief Get min error of the linear solver for the poisson equation. + * \return Min error of the linear solver for the poisson equation. + */ + su2double GetPoisson_Linear_Solver_Error(void) const { return Poisson_Linear_Solver_Error; } + /*! * \brief Get max number of iterations of the linear solver for the implicit formulation. * \return Max number of iterations of the linear solver for the implicit formulation. */ unsigned long GetLinear_Solver_Iter(void) const { return Linear_Solver_Iter; } + /*! + * \brief Get max number of iterations of the linear solver for the poisson equation. + * \return Max number of iterations of the linear solver for the poisson equation. + */ + unsigned long GetPoisson_Linear_Solver_Iter(void) const { return Poisson_Linear_Solver_Iter; } + /*! * \brief Get max number of iterations of the linear solver for the implicit formulation. * \return Max number of iterations of the linear solver for the implicit formulation. @@ -4392,6 +4444,11 @@ class CConfig { */ const CIluOptions& GetIluOptions(void) const { return IluOptions; } + /*! + * \brief Get the SIMPLE (and PISO) algorithm options, see CSIMPLE_Options. + */ + const CSIMPLE_Options& GetSIMPLE_Options(void) const { return SIMPLE_Options; } + /*! * \brief Get restart frequency of the linear solver for the implicit formulation. * \return Restart frequency of the linear solver for the implicit formulation. diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index abecc767e26..ca4be5e1a91 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -38,6 +38,7 @@ #include #include #include +#include /*--- In forward mode the matrix is not of a built-in type. ---*/ #if defined(HAVE_MKL) && !defined(CODI_FORWARD_TYPE) @@ -678,15 +679,15 @@ class CSysMatrix { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \param[in] needTranspPtr - If the L/U transpose maps should be built, used for "SetDiagonalAsColumnSum". - * \param[in] grad_mode - Gradient smoothing mode, only used to detect the right preconditioner type. * \param[in] allow_quant - Quantization is only possible with solvers that "set and forget" the off-diagonal * blocks of the matrix. Solvers that perform multiple updates would lose too much information, so * that pattern is not supported with quantization (the code will hit null pointers). It is up to * the solver to declare whether it will "set and forget". + * \param[in] override_prec - Decide if, and with what argument to override the preconditioner. */ void Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry* geometry, const CConfig* config, bool needTranspPtr = false, - bool grad_mode = false, bool allow_quant = false); + bool allow_quant = false, std::optional override_prec = std::nullopt); /*! * \brief Compresses off-diagonal blocks into quantized form for use with USE_QUANTIZATION. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index d4d0e7f0fa0..0ce98fd04b7 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -77,7 +77,7 @@ const unsigned int MAX_PARAMETERS = 10; /*!< \brief Maximum number of para const unsigned int MAX_NUMBER_PERIODIC = 10; /*!< \brief Maximum number of periodic boundary conditions. */ const unsigned int MAX_STRING_SIZE = 400; /*!< \brief Maximum size of a generic string. */ const unsigned int MAX_NUMBER_FFD = 15; /*!< \brief Maximum number of FFDBoxes for the FFD. */ -enum: unsigned int{MAX_SOLS = 13}; /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */ +enum: unsigned int{MAX_SOLS = 14}; /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */ const unsigned int MAX_TERMS = 7; /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */ const unsigned int MAX_ZONES = 3; /*!< \brief Maximum number of zones. */ const unsigned int MAX_FE_KINDS = 7; /*!< \brief Maximum number of Finite Elements. */ @@ -264,6 +264,7 @@ enum class MAIN_SOLVER { FEM_RANS, /*!< \brief Definition of the finite element Reynolds-averaged Navier-Stokes' (RANS) solver. */ FEM_LES, /*!< \brief Definition of the finite element Large Eddy Simulation Navier-Stokes' (LES) solver. */ MULTIPHYSICS, + POISSON_EQUATION, /*!< \brief Definition of the Poisson equation solver. */ NEMO_EULER, /*!< \brief Definition of the NEMO Euler solver. */ NEMO_NAVIER_STOKES, /*!< \brief Definition of the NEMO NS solver. */ }; @@ -339,6 +340,31 @@ static const MapType MatComp_Map = { MakePair("NEARLY_INCOMPRESSIBLE", STRUCT_COMPRESS::NEARLY_INCOMP) }; +/*! + * \brief Type of incompressible solver + */ +enum class INCOMP_SYSTEM { + DENSITY_BASED, /*!< \brief Density-based. */ + PRESSURE_BASED, /*!< \brief Pressure-based. */ +}; +static const MapType Incomp_Map = { + MakePair("DENSITY_BASED", INCOMP_SYSTEM::DENSITY_BASED) + MakePair("PRESSURE_BASED", INCOMP_SYSTEM::PRESSURE_BASED) +}; + +/*! + * \brief Type of iteration + */ +enum class PBITER { + SIMPLE, /*!< \brief SIMPLE algorithm. */ + SIMPLEC, /*!< \brief SIMPLEC algorithm. */ +}; + +static const MapType PBIter_Map = { + MakePair("SIMPLE", PBITER::SIMPLE) + MakePair("SIMPLEC", PBITER::SIMPLEC) +}; + /*! * \brief Types of interpolators */ @@ -479,6 +505,7 @@ enum RUNTIME_TYPE { RUNTIME_ADJRAD_SYS = 24, /*!< \brief One-physics case, the code is solving the adjoint radiation model. */ RUNTIME_SPECIES_SYS = 25, /*!< \brief One-physics case, the code is solving the species model. */ RUNTIME_ADJSPECIES_SYS = 26,/*!< \brief One-physics case, the code is solving the adjoint species model. */ + RUNTIME_POISSON_SYS = 27, /*!< \brief One-physics case, the code is solving the poisson equation. */ }; enum SOLVER_TYPE : const int { @@ -497,6 +524,7 @@ enum RUNTIME_TYPE { ADJSPECIES_SOL=12, /*!< \brief Position of the adjoint of the species solver. */ FEA_SOL=0, /*!< \brief Position of the Finite Element flow solution in the solver container array. */ ADJFEA_SOL=1, /*!< \brief Position of the continuous adjoint Finite Element flow solution in the solver container array. */ + POISSON_SOL=13, /*!< \brief Position of the poisson solution in the solver container array */ TEMPLATE_SOL=0, /*!< \brief Position of the template solution. */ }; @@ -828,7 +856,8 @@ enum class CENTERED { LAX, /*!< \brief Lax-Friedrich centered numerical method. */ JST_MAT, /*!< \brief JST with matrix dissipation. */ JST_KE, /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */ - LD2 /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */ + LD2, /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */ + CDS /*!< \brief Central Difference Scheme used for pressure based solver. */ }; static const MapType Centered_Map = { MakePair("NONE", CENTERED::NONE) @@ -837,6 +866,7 @@ static const MapType Centered_Map = { MakePair("JST_MAT", CENTERED::JST_MAT) MakePair("LAX-FRIEDRICH", CENTERED::LAX) MakePair("LD2", CENTERED::LD2) + MakePair("CDS", CENTERED::CDS) }; @@ -863,7 +893,8 @@ enum class UPWIND { AUSMPLUSUP, /*!< \brief AUSM+ -up numerical method (All Speed) */ AUSMPLUSUP2, /*!< \brief AUSM+ -up2 numerical method (All Speed) */ AUSMPLUSM, /*!< \breif AUSM+M numerical method. (NEMO Only)*/ - BOUNDED_SCALAR /*!< \brief Scalar advection numerical method. */ + BOUNDED_SCALAR, /*!< \brief Scalar advection numerical method. */ + UDS /*!< \brief Upwind Difference Scheme used for pressure based solver. */ }; static const MapType Upwind_Map = { MakePair("NONE", UPWIND::NONE) @@ -885,6 +916,7 @@ static const MapType Upwind_Map = { MakePair("SLAU2", UPWIND::SLAU2) MakePair("FDS", UPWIND::FDS) MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH) + MakePair("UDS", UPWIND::UDS) }; /*! @@ -2781,6 +2813,9 @@ enum class MPI_QUANTITIES { MESH_DISPLACEMENTS , /*!< \brief Mesh displacements at the interface. */ SOLUTION_TIME_N , /*!< \brief Solution at time n. */ SOLUTION_TIME_N1 , /*!< \brief Solution at time n-1. */ + MOM_COEFF , /*!< \brief Momentum coefficient for the Rhie-Chow scheme. */ + MOM_CORRECTION , /*!< \brief Momentum correction for the pressure-based poisson solver (used when computing HbyA). */ + HBYA_CORRECTION , /*!< \brief HbyA correction for the pressure-based poisson solver. */ }; /*! @@ -2908,6 +2943,7 @@ enum class LINEAR_SOLVER_MODE { STANDARD, /*!< \brief Operate in standard mode. */ MESH_DEFORM, /*!< \brief Operate in mesh deformation mode. */ GRADIENT_MODE, /*!< \brief Operate in gradient smoothing mode. */ + POISSON, /*!< \brief Operate in poisson solver mode. */ }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c21116c6f3e..2a0b9d3d021 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1202,6 +1202,11 @@ void CConfig::SetConfig_Options() { /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map); + /*!\brief KIND_INCOMP_SYSTEM \n DESCRIPTION: Incomp type \n OPTIONS: see \link Incomp_Map \endlink DEFAULT: NONE \ingroup Config*/ + addEnumOption("KIND_INCOMP_SYSTEM", Kind_Incomp_System, Incomp_Map, INCOMP_SYSTEM::DENSITY_BASED); + /*!\brief KIND_PB_ITER \n DESCRIPTION: Kind_PBIter \n OPTIONS: see \link PBIter_Map \endlink \ingroup Config*/ + addEnumOption("KIND_PB_ITER", Kind_PBIter, PBIter_Map, PBITER::SIMPLE); + /*!\brief ROUGHSST_OPTIONS \n DESCRIPTION: Specify type of boundary condition for rough walls for SST turbulence model. \n Options: see \link ROUGHSST_Options_Map \endlink \n DEFAULT: wilcox1998 \ingroup Config*/ addEnumOption("KIND_ROUGHSST_MODEL", Kind_RoughSST_Model, RoughSST_Model_Map, ROUGHSST_MODEL::WILCOX1998); /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ @@ -1968,6 +1973,16 @@ void CConfig::SetConfig_Options() { addDoubleOption("LINEAR_SOLVER_ERROR", Linear_Solver_Error, 1E-6); /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */ addUnsignedLongOption("LINEAR_SOLVER_ITER", Linear_Solver_Iter, 10); + /*!\brief LINEAR_SOLVER + * \n DESCRIPTION: Linear solver for the poisson system \n OPTIONS: see \link Linear_Solver_Map \endlink \n DEFAULT: FGMRES \ingroup Config*/ + addEnumOption("POISSON_LINEAR_SOLVER", Kind_Poisson_Linear_Solver, Linear_Solver_Map, FGMRES); + /*!\brief LINEAR_SOLVER_PREC + * \n DESCRIPTION: Preconditioner for the Krylov linear solvers \n OPTIONS: see \link Linear_Solver_Prec_Map \endlink \n DEFAULT: LU_SGS \ingroup Config*/ + addEnumOption("POISSON_LINEAR_SOLVER_PREC", Kind_Poisson_Linear_Solver_Prec, Linear_Solver_Prec_Map, ILU); + /* DESCRIPTION: Minimum error threshold for the poisson linear solver */ + addDoubleOption("POISSON_LINEAR_SOLVER_ERROR", Poisson_Linear_Solver_Error, 1E-6); + /* DESCRIPTION: Maximum number of iterations of the poisson linear solver */ + addUnsignedLongOption("POISSON_LINEAR_SOLVER_ITER", Poisson_Linear_Solver_Iter, 10); /* DESCRIPTION: Fill in level for the ILU preconditioner */ addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", IluOptions.FillIn, 0); /* DESCRIPTION: Use level scheduling for OMP parallelization of the ILU preconditioner */ @@ -1984,6 +1999,14 @@ void CConfig::SetConfig_Options() { addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0); /* DESCRIPTION: Use an inner linear solver. */ addEnumOption("LINEAR_SOLVER_INNER", Kind_Linear_Solver_Inner, Inner_Linear_Solver_Map, LINEAR_SOLVER_INNER::NONE); + /* DESCRIPTION: Relaxation of the pressure corrections for the SIMPLE algorithm */ + addDoubleOption("RELAXATION_FACTOR_PRESSURE", SIMPLE_Options.Relaxation_Factor_Pressure, 1.0); + /* DESCRIPTION: Removal factor for the transient term in the momentum coefficients for the poisson solver. */ + addDoubleOption("TRANSIENT_TERM_REMOVAL_FACTOR", SIMPLE_Options.Transient_Term_Removal_Factor, 0.0); + /*!\DESCRIPTION: Automatically compute relaxation factors for flow corrections in the SIMPLE algorithm */ + addBoolOption("USE_AUTOMATIC_RELAXATION_FACTORS", SIMPLE_Options.AutomaticRelaxationFactors, false); + /* DESCRIPTION: Number of corrections in the PISO algorithm (pressure based). */ + addUnsignedShortOption("PISO_CORRECTIONS", SIMPLE_Options.nCorrections_PISO, 1); /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */ addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0); /* DESCRIPTION: Relaxation of the CHT coupling */ @@ -8939,6 +8962,7 @@ unsigned short CConfig::GetContainerPosition(unsigned short val_eqsystem) { case RUNTIME_ADJSPECIES_SYS:return ADJSPECIES_SOL; case RUNTIME_ADJFEA_SYS: return ADJFEA_SOL; case RUNTIME_RADIATION_SYS: return RAD_SOL; + case RUNTIME_POISSON_SYS: return POISSON_SOL; case RUNTIME_MULTIGRID_SYS: return 0; } return 0; @@ -9077,6 +9101,13 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, } break; + case MAIN_SOLVER::POISSON_EQUATION: + if (val_system == RUNTIME_POISSON_SYS) { + SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, 0.0, NONE); + SetKind_TimeIntScheme(EULER_IMPLICIT); + } + break; + case MAIN_SOLVER::FEM_ELASTICITY: case MAIN_SOLVER::DISC_ADJ_FEM: if (val_system == RUNTIME_FEA_SYS) { diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 84f30a2aa3d..7c26da965b9 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -179,7 +179,8 @@ CSysMatrix::~CSysMatrix() { template void CSysMatrix::Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry* geometry, - const CConfig* config, bool needTranspPtr, bool grad_mode, bool allow_quant) { + const CConfig* config, bool needTranspPtr, bool allow_quant, + std::optional override_prec) { SU2_ZONE_SCOPED assert(omp_get_thread_num() == 0 && "Only the master thread is allowed to initialize the matrix."); @@ -209,8 +210,8 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi } /*--- No else if, but separate if case! ---*/ - if (config->GetSmoothGradient() && grad_mode) { - prec = config->GetKind_Grad_Linear_Solver_Prec(); + if (override_prec) { + prec = *override_prec; } useCuda = config->GetCUDA(); diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 3c9176d8716..cd88e8ce5ce 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -1421,6 +1421,16 @@ unsigned long CSysSolve::Solve(CSysMatrix& Jacobian, con break; } + /*--- Poisson solver mode ---*/ + case LINEAR_SOLVER_MODE::POISSON: { + KindSolver = config->GetKind_Poisson_Linear_Solver(); + KindPrecond = config->GetKind_Poisson_Linear_Solver_Prec(); + MaxIter = config->GetPoisson_Linear_Solver_Iter(); + SolverTol = SU2_TYPE::GetValue(config->GetPoisson_Linear_Solver_Error()); + ScreenOutput = false; + break; + } + /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', * but does not enforce it to avoid compiler warning. ---*/ default: { @@ -1632,6 +1642,16 @@ unsigned long CSysSolve::Solve_b(CSysMatrix& Jacobian, c break; } + /*--- Poisson solver mode ---*/ + case LINEAR_SOLVER_MODE::POISSON: { + KindSolver = config->GetKind_Poisson_Linear_Solver(); + KindPrecond = config->GetKind_Poisson_Linear_Solver_Prec(); + MaxIter = config->GetPoisson_Linear_Solver_Iter(); + SolverTol = SU2_TYPE::GetValue(config->GetPoisson_Linear_Solver_Error()); + ScreenOutput = false; + break; + } + /*--- Normal mode assumes that 'lin_sol_mode==LINEAR_SOLVER_MODE::STANDARD', * but does not enforce it to avoid compiler warning. ---*/ default: { diff --git a/SU2_CFD/include/iteration/CFluidIteration.hpp b/SU2_CFD/include/iteration/CFluidIteration.hpp index 024f88e7d43..ec0ca704cb6 100644 --- a/SU2_CFD/include/iteration/CFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CFluidIteration.hpp @@ -72,6 +72,27 @@ class CFluidIteration : public CIteration { CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) override; + /*! + * \brief Perform auxiliary solvers iterations after the main flow solver. + * \param[in] output - Pointer to the COutput class. + * \param[in] integration - Container vector with all the integration methods. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method (the way in which the equations are solved). + * \param[in] config - Definition of the particular problem. + * \param[in] surface_movement - Surface movement classes of the problem. + * \param[in] grid_movement - Volume grid movement classes of the problem. + * \param[in] FFDBox - FFD FFDBoxes of the problem. + * \param[in] val_iZone - Index of the zone. + * \param[in] val_iInst - Index of the instance layer. + * \param[in] main_solver - Main solver. + * \param[in] frozen_visc - Flag for frozen viscosity. + */ + void CommonAuxiliarySolvers(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, + CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, + CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, + unsigned short val_iInst, MAIN_SOLVER main_solver, bool frozen_visc); + /*! * \brief Iterate the fluid system for a number of Inner_Iter iterations. * \param[in] output - Pointer to the COutput class. diff --git a/SU2_CFD/include/iteration/CPBFluidIteration.hpp b/SU2_CFD/include/iteration/CPBFluidIteration.hpp new file mode 100644 index 00000000000..92f0c16a49f --- /dev/null +++ b/SU2_CFD/include/iteration/CPBFluidIteration.hpp @@ -0,0 +1,66 @@ +/*! + * \file CPBFluidIteration.hpp + * \brief Headers of the pressure based fluid iteration class. + * \author F. Palacios, T. Economon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CFluidIteration.hpp" + +/*! + * \class CFluidIteration + * \ingroup Drivers + * \brief Class for driving an iteration of the fluid system. + * \author T. Economon + */ +class CPBFluidIteration : public CFluidIteration { +public: + /*! + * \brief Constructor of the class. + * \param[in] config - Definition of the particular problem. + */ + explicit CPBFluidIteration(const CConfig* config) : CFluidIteration(config) {} + + /*! + * \brief Perform a single iteration of the fluid system. + * \param[in] output - Pointer to the COutput class. + * \param[in] integration - Container vector with all the integration methods. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method (the way in which the equations are solved). + * \param[in] config - Definition of the particular problem. + * \param[in] surface_movement - Surface movement classes of the problem. + * \param[in] grid_movement - Volume grid movement classes of the problem. + * \param[in] FFDBox - FFD FFDBoxes of the problem. + * \param[in] val_iZone - Index of the zone. + * \param[in] val_iInst - Index of the instance layer. + */ + void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, + CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, + CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, + unsigned short val_iInst) override; + + }; + diff --git a/SU2_CFD/include/numerics/flow/convection/pressure_based.hpp b/SU2_CFD/include/numerics/flow/convection/pressure_based.hpp new file mode 100644 index 00000000000..6e1d67303c3 --- /dev/null +++ b/SU2_CFD/include/numerics/flow/convection/pressure_based.hpp @@ -0,0 +1,160 @@ +/*! + * \file pressure_based.hpp + * \brief Declaration of numerics classes for convective schemes for + * the pressure based solver, the implementation is in pressure_based.cpp. + * \author T. Aalbers + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../../CNumerics.hpp" + + +/*! + * \class CPBConvection_Base + * \brief Class for computing a linear centered scheme. + * \ingroup ConvDiscr + * \author T. Aalbers + */ +class CPBConvection_Base : public CNumerics { +protected: + + bool implicit, dynamic_grid, energy, variable_density; + + unsigned short iDim, jDim, iVar, jVar; + + su2double *AdvectedVelocity = nullptr, AdvectedEnthalpy; + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + + su2double MeanPressure, MeanDensity; + su2double dRhodh_i, dRhodh_j, Temperature_i, Temperature_j; + + su2double weight_jacobian_i, weight_jacobian_j; + + /*! + * \brief Function which defines the advected quantities + */ + void virtual ComputeAdvectedQuantities(void) = 0; + + /*! + * \brief Function which defines jacobian weights + */ + void virtual ComputeJacobianWeights(void) = 0; + + /*! + * \brief Function which defines the Jacobian + */ + void ComputeJacobian(su2double val_density, const su2double *val_velocity, + su2double val_enthalpy, su2double val_dRhodh, + su2double val_scale, su2double **val_Proj_Jac_Tensor); + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimension of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CPBConvection_Base(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + virtual ~CPBConvection_Base(void); + + /*! + * \brief Compute the flow residual. + * \param[out] val_resconv - Pointer to the convective residual. + * \param[out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[out] val_Jacobian_j - Jacobian of the numerical method at node j (implicit computation). + * \param[in] config - Definition of the particular problem. + */ + ResidualType<> ComputeResidual(const CConfig* config) final; +}; + +/*! + * \class CPBConvection_Central + * \brief Class for computing a centered scheme. + * \ingroup ConvDiscr + * \author T. Aalbers + */ +class CPBConvection_Central : public CPBConvection_Base { + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimension of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CPBConvection_Central(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) + : CPBConvection_Base(val_nDim, val_nVar, config) {} + + /*! + * \brief Function which defines the advected quantities + */ + void ComputeAdvectedQuantities(void) final; + + /*! + * \brief Function which defines jacobian weights + */ + void ComputeJacobianWeights(void) final; + +}; + + +/*! + * \class CPBConvection_Upwind + * \brief Class for computing an upwind scheme. + * \ingroup ConvDiscr + * \author T. Aalbers + */ +class CPBConvection_Upwind : public CPBConvection_Base { + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimension of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CPBConvection_Upwind(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) + : CPBConvection_Base(val_nDim, val_nVar, config) {} + + /*! + * \brief Function which defines the advected quantities + */ + void ComputeAdvectedQuantities(void) final; + + /*! + * \brief Function which defines jacobian weights + */ + void ComputeJacobianWeights(void) final; + +}; \ No newline at end of file diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp index 52c1d6e743e..f5ce83935a1 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp @@ -82,8 +82,8 @@ CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CV case CENTERED::JST_MAT: obj = new CJSTmatScheme(config, iMesh, turbVars); break; - case CENTERED::LD2: - /*--- LD2 implemented only in the incompressible solver. ---*/ + case CENTERED::LD2: case CENTERED::CDS: + /*--- CDS and LD2 implemented only in the incompressible solver. ---*/ break; } return obj; diff --git a/SU2_CFD/include/output/CFlowIncOutput.hpp b/SU2_CFD/include/output/CFlowIncOutput.hpp index 0eb29f2911a..80e22a80d3b 100644 --- a/SU2_CFD/include/output/CFlowIncOutput.hpp +++ b/SU2_CFD/include/output/CFlowIncOutput.hpp @@ -42,6 +42,7 @@ class CFlowIncOutput final: public CFlowOutput { bool heat; /*!< \brief Boolean indicating whether have a heat problem*/ bool weakly_coupled_heat; /*!< \brief Boolean indicating whether have a weakly coupled heat equation*/ bool flamelet; /*!< \brief Boolean indicating whether we solve the flamelet equations */ + bool pressure_based; /*!< Boolean indicating whether running the pressure based version */ unsigned short streamwisePeriodic; /*!< \brief Boolean indicating whether it is a streamwise periodic simulation. */ bool streamwisePeriodic_temperature; /*!< \brief Boolean indicating streamwise periodic temperature is used. */ diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 9c94d4195e9..fd5fbf98609 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -41,6 +41,12 @@ class CIncEulerSolver : public CFVMFlowSolverBase FluidModel; /*!< \brief fluid model used in the solver. */ StreamwisePeriodicValues SPvals, SPvalsUpdated; + bool pressure_based; + su2activevector alpha_p; + su2activevector pressureCorrection; + su2activematrix momentumCorrection; + su2activevector EdgeMassFluxCorrection; + /*! * \brief Preprocessing actions common to the Euler and NS solvers. * \param[in] geometry - Geometrical definition of the problem. @@ -148,6 +154,22 @@ class CIncEulerSolver : public CFVMFlowSolverBase. + */ + +#pragma once + +#include "CScalarSolver.hpp" +#include "../variables/CPoissonVariable.hpp" + +/*! + * \class CPoissonSolver + * \brief Main class for defining the finite-volume poisson equation solver. + * \author O. Burghardt + * \version 8.5.0 "Harrier" + */ +class CPoissonSolver final : public CScalarSolver { +protected: + static constexpr size_t MAXNDIM = 3; /*!< \brief Max number of space dimensions, used in some static arrays. */ + static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */ + + /*! + * \brief Compute the viscous flux for the scalar equation at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \note Calls a generic implementation after defining a SolverSpecificNumerics object. + */ + inline void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, + CNumerics* numerics, const CConfig* config) override { + + su2double mom_coeff_i{}, mom_coeff_j{}; + + /*--- Sets the momentum coefficients to use in the viscous numerics. ---*/ + auto compute_momentum_coeff = [&](unsigned long iPoint, unsigned long jPoint) { + mom_coeff_i = nodes->GetMomCoeff(iPoint); + mom_coeff_j = nodes->GetMomCoeff(jPoint); + numerics->SetDiffusionCoeff(&mom_coeff_i, &mom_coeff_j); + }; + + /*--- Compute residual and Jacobians. ---*/ + Viscous_Residual_impl(compute_momentum_coeff, iEdge, geometry, solver_container, numerics, config); + } + +public: + + /* + * \overload + * \param[in] geometry - Geometrical definition of the problem + * \param[in] config - Definition of the particular problem + */ + CPoissonSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh); + + /*! + * \brief Restart residual and compute gradients. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void Preprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned short iRKStep, + unsigned short RunTime_EqSystem, + bool Output) override; + + /*! + * \brief Correct the pressure and velocities for the flow solution + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Postprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh) final; + + /*! + * \brief Compute the viscous residuals for the turbulent equation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics_container - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + */ + void Viscous_Residual(CGeometry *geometry, + CSolver **solver_container, + CNumerics **numerics_container, + CConfig *config, + unsigned short iMesh, + unsigned short iRKStep) override; + + /*! + * \brief Source term computation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics_container - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, + CConfig *config, unsigned short iMesh) override; + + /*! + * \brief No upwind residual for poisson equation. + */ + void Upwind_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, + CConfig* config, unsigned short iMesh) override {} + + /*! + * \brief Update the solution using an implicit solver. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + void ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) override; + + /*! + * \brief No time step as it is a stationary problem + */ + void SetTime_Step(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned long Iteration) override {} + + /*! + * \brief No dual time stepping as there is no time stepping at all. + */ + void SetResidual_DualTime(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iRKStep, + unsigned short iMesh, unsigned short RunTime_EqSystem) override {} + + /*! + * \brief Compute the coefficients for the pressure correction equation based + * on the residuals from the solution of the momentum equation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] periodic - Flag for periodic boundary conditions. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void SetMomCoeff(CGeometry *geometry, CSolver **solver_container, CConfig *config, bool periodic, unsigned short iMesh) final; + + + /*! + * \brief Compute the HbyA based on the momentum correction to be used in second PISO + * correction equation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void ComputeHbyA(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) final; + + /*! + * \brief Impose a constant heat-flux condition at the wall. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_HeatFlux_Wall(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) final; + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, + unsigned short val_marker) final; + + /*! + * \brief Impose the inlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Inlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + /*! + * \brief Impose the outlet boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Outlet(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + +}; diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index 500233645bf..b09ff29ee61 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -441,7 +441,7 @@ class CScalarSolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar); + CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar, LINEAR_SOLVER_MODE linear_solver_mode = LINEAR_SOLVER_MODE::STANDARD); /*! * \brief Compute the spatial integration using a upwind scheme. @@ -464,7 +464,7 @@ class CScalarSolver : public CSolver { * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) final; + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) override; /*! * \brief Impose the Symmetry Plane boundary condition. @@ -602,7 +602,7 @@ class CScalarSolver : public CSolver { * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) final; + void ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) override; /*! * \brief Set the total residual adding the term that comes from the Dual Time-Stepping Strategy. @@ -625,7 +625,7 @@ class CScalarSolver : public CSolver { * \param[in] val_update_geo - Flag for updating coords and grid velocity. */ void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, - bool val_update_geo) override = 0; + bool val_update_geo) override {} /*! * \brief Scalar solvers support OpenMP+MPI. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 2a20484449b..a15d068bbb5 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -30,8 +30,8 @@ #include "../../include/variables/CFlowVariable.hpp" template -CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar) - : CSolver(), Conservative(conservative), BoundedScalar(bounded_scalar), +CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar, LINEAR_SOLVER_MODE linear_solver_mode) + : CSolver(linear_solver_mode), Conservative(conservative), BoundedScalar(bounded_scalar), prim_idx(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE, config->GetNEMOProblem(), geometry->GetnDim(), config->GetnSpecies()) { SU2_ZONE_SCOPED diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 863eb77b49f..332eba9a4bd 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -777,6 +777,22 @@ class CSolver { unsigned short RunTime_EqSystem, bool Output) { } + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + inline virtual void ComputeEdgeMassFluxesRhieChow(CGeometry *geometry, CSolver **solver_container, CConfig *config) { } + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + inline virtual void ApplyPressureVelocityCorrection(CGeometry *geometry, CSolver **solver_container, CConfig *config) { } + /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -4298,6 +4314,25 @@ class CSolver { */ virtual StreamwisePeriodicValues GetStreamwisePeriodicValues() const { return StreamwisePeriodicValues(); } + /*! + * \brief A virtual member + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] periodic - Flag for periodic boundary conditions. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + inline virtual void SetMomCoeff(CGeometry *geometry, CSolver **solver_container, CConfig *config, bool periodic, unsigned short iMesh) { } + + /*! + * \brief A virtual member + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + inline virtual void ComputeHbyA(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { } + /*! * \brief Save snapshot or POD data using libROM * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CSolverFactory.hpp b/SU2_CFD/include/solvers/CSolverFactory.hpp index e612f3fb1e3..9a89629658e 100644 --- a/SU2_CFD/include/solvers/CSolverFactory.hpp +++ b/SU2_CFD/include/solvers/CSolverFactory.hpp @@ -63,6 +63,7 @@ enum class SUB_SOLVER_TYPE { MESH, /*!< \brief Mesh solver */ RADIATION, /*!< \brief Radiation solver */ DISC_ADJ_RADIATION, /*!< \brief Discrete adjoint radiation solver */ + POISSON, /*!< \breif Poisson equation solver */ NONE }; diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index d76990b71eb..f039a112779 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -74,6 +74,10 @@ class CIncEulerVariable : public CFlowVariable { VectorType Density_time_n, /*!< \brief Density at time n for dual-time stepping. */ Density_time_n1; /*!< \brief Density at time n-1 for dual-time stepping. */ su2double TemperatureLimits[2]; /*!< \brief Temperature limits [K]. */ + + using BoolVectorType = C2DContainer; + BoolVectorType strongBC; /*!< \brief Flag for boundary conditions to indicate if a strong BC has been applied, currently only used to keep track of farfield. */ + public: /*! * \brief Constructor of the class. @@ -325,4 +329,23 @@ class CIncEulerVariable : public CFlowVariable { */ inline void SetDensity_time_n1(unsigned long iPoint, su2double val_density) { Density_time_n1(iPoint) = val_density; } + /*! + * \brief Set the BC flag to true of the point. + * \param[in] iPoint - Point index. + */ + inline void SetStrongBC(unsigned long iPoint) { strongBC(iPoint) = true; } + + /*! + * \brief Get the BC flag of the point + * \param[in] iPoint - Point index. + * \return The boolean flag of the strong boundary condition. + */ + inline bool GetStrongBC(unsigned long iPoint) { return strongBC(iPoint); } + + /*! + * \brief Set the BC flag to false of the point. + * \param[in] iPoint - Point index. + */ + inline void ResetStrongBC(unsigned long iPoint) { strongBC(iPoint) = false; } + }; diff --git a/SU2_CFD/include/variables/CPoissonVariable.hpp b/SU2_CFD/include/variables/CPoissonVariable.hpp new file mode 100644 index 00000000000..6dfdc839adc --- /dev/null +++ b/SU2_CFD/include/variables/CPoissonVariable.hpp @@ -0,0 +1,99 @@ +/*! + * \file CPoissonVariable.hpp + * \brief Class for defining the variables of the finite-volume heat equation solver. + * \author F. Palacios, T. Economon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "CScalarVariable.hpp" + +/*! + * \class CPoissonVariable + * \brief Class for defining the variables of the finite-volume poisson equation solver. + * \author O. Burghardt + * \version 8.5.0 "Harrier" + */ +class CPoissonVariable final : public CScalarVariable { +protected: + VectorType MomCoeff; /*!< \brief Momentum coefficients vol/A_p used as the diffusion coefficients in the poisson solver. */ + MatrixType MomentumCorrection; /*!< \brief (rho*u)' in the context of: (rho*u)** = (rho*u)* + (rho*u)'. */ + MatrixType HbyACorrection; /*!< \brief H(rhou')/A = (sum_nb A_nb (rhou)'_nb) / A; used by the second pressure correction in the PISO algorithm. */ +public: + static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */ + + /*! + * \brief Constructor of the class. + * \param[in] value - Values of the poisson solution (initialization value). + * \param[in] npoint - Number of points/nodes/vertices in the domain. + * \param[in] ndim - Number of dimensions of the problem. + * \param[in] nvar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CPoissonVariable(su2double value, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config); + + /*! + * \brief Get the momentum coefficient of the point. + * \return Value of the momentum coefficient of the point. + */ + inline su2double GetMomCoeff(unsigned long iPoint) final { return MomCoeff(iPoint);} + + /*! + * \brief Set the momentum coefficient of the point. + */ + inline void SetMomCoeff(unsigned long iPoint, su2double val_Mom_Coeff) final { MomCoeff(iPoint) = val_Mom_Coeff; } + + /*! + * \brief Set H(u')/A for the point + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_HbyA - HbyA correction. + */ + inline void SetHbyACorrection(unsigned long iPoint, unsigned short iDim, su2double val_HbyA) final { HbyACorrection(iPoint, iDim) = val_HbyA; } + + /*! + * \brief Get H(u')/A for the point + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \return The H(u')/A for the point. + */ + inline su2double GetHbyACorrection(unsigned long iPoint, unsigned short iDim) final { return HbyACorrection(iPoint, iDim); } + + /*! + * \brief Set (rho*u)' for the point + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_mom - Momentum correction (rho*u)' value. + */ + inline void SetMomentumCorrection(unsigned long iPoint, unsigned short iDim, su2double val_mom) final { MomentumCorrection(iPoint, iDim) = val_mom; } + + /*! + * \brief Get (rho*u)' for the point + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \return The (rho*u)' for the point. + */ + inline su2double GetMomentumCorrection(unsigned long iPoint, unsigned short iDim) final { return MomentumCorrection(iPoint, iDim); } + +}; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 528adc139db..21d50420380 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2416,4 +2416,16 @@ class CVariable { inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; } inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; } + + inline virtual su2double GetMomCoeff(unsigned long iPoint) { return 0.0; } + + inline virtual void SetMomCoeff(unsigned long iPoint, su2double val_Mom_Coeff) { } + + inline virtual su2double GetMomentumCorrection(unsigned long iPoint, unsigned short iDim) { return 0.0; } + + inline virtual void SetMomentumCorrection(unsigned long iPoint, unsigned short iDim, su2double val_mom) { } + + inline virtual su2double GetHbyACorrection(unsigned long iPoint, unsigned short iDim) { return 0.0; } + + inline virtual void SetHbyACorrection(unsigned long iPoint, unsigned short iDim, su2double val_HbyAcorrection) { } }; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 5904ed9febb..3f53610dc35 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -62,6 +62,7 @@ #include "../../include/numerics/flow/convection/hllc.hpp" #include "../../include/numerics/flow/convection/ausm_slau.hpp" #include "../../include/numerics/flow/convection/centered.hpp" +#include "../../include/numerics/flow/convection/pressure_based.hpp" #include "../../include/numerics/flow/flow_diffusion.hpp" #include "../../include/numerics/flow/flow_sources.hpp" #include "../../include/numerics/NEMO/convection/roe.hpp" @@ -1455,6 +1456,8 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver bool compressible = false; bool incompressible = false; bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS); + bool pressure_based = (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED); + bool poisson = (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED); bool roe_low_dissipation = (config->GetKind_RoeLowDiss() != NO_ROELOWDISS); /*--- Initialize some useful booleans ---*/ @@ -1567,7 +1570,7 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver if (fem_ns) nVar_Flow = solver[MESH_0][FLOW_SOL]->GetnVar(); if (fem) nVar_FEM = solver[MESH_0][FEA_SOL]->GetnVar(); - + if (config->AddRadiation()) nVar_Rad = solver[MESH_0][RAD_SOL]->GetnVar(); /*--- Number of variables for adjoint problem ---*/ @@ -1655,22 +1658,36 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } if (incompressible) { - /*--- Incompressible flow, use preconditioning method ---*/ - switch (config->GetKind_Centered_Flow()) { - case CENTERED::LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break; - case CENTERED::LD2 : - case CENTERED::JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break; - default: - SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only JST and LAX-FRIEDRICH are available for incompressible flows.", CURRENT_FUNCTION); - break; - } - for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); + if (!pressure_based) { + /*--- Incompressible flow, use preconditioning method ---*/ + switch (config->GetKind_Centered_Flow()) { + case CENTERED::LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break; + case CENTERED::LD2 : + case CENTERED::JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break; + default: + SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only JST and LAX-FRIEDRICH are available for density based incompressible flows.", CURRENT_FUNCTION); + break; + } + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); + /*--- Definition of the boundary condition method ---*/ + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); - /*--- Definition of the boundary condition method ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); + } else { + /*--- Incompressible flow, use pressure-based method ---*/ + switch (config->GetKind_Centered_Flow()) { + case CENTERED::CDS : numerics[MESH_0][FLOW_SOL][conv_term] = new CPBConvection_Central(nDim, nVar_Flow, config); break; + default: + SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only CDS is available for pressure based incompressible flows.", CURRENT_FUNCTION); + } + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][conv_term] = new CPBConvection_Central(nDim, nVar_Flow, config); + /*--- Definition of the boundary condition method ---*/ + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CPBConvection_Upwind(nDim, nVar_Flow, config); + } } break; case SPACE_UPWIND : @@ -1777,17 +1794,32 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } if (incompressible) { - /*--- Incompressible flow, use artificial compressibility method ---*/ - switch (config->GetKind_Upwind_Flow()) { - case UPWIND::FDS: - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); - numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); - } - break; - default: - SU2_MPI::Error("Invalid upwind scheme or not implemented.\n Currently, only FDS is available for incompressible flows.", CURRENT_FUNCTION); - break; + if (!pressure_based) { + /*--- Incompressible flow, use artificial compressibility method ---*/ + switch (config->GetKind_Upwind_Flow()) { + case UPWIND::FDS: + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); + numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwFDSInc_Flow(nDim, nVar_Flow, config); + } + break; + default: + SU2_MPI::Error("Invalid upwind scheme or not implemented.\n Currently, only FDS is available for density based incompressible flows.", CURRENT_FUNCTION); + break; + } + } else { + /*--- Incompressible flow, use pressure based method ---*/ + switch (config->GetKind_Upwind_Flow()) { + case UPWIND::UDS: + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + numerics[iMGlevel][FLOW_SOL][conv_term] = new CPBConvection_Upwind(nDim, nVar_Flow, config); + numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CPBConvection_Upwind(nDim, nVar_Flow, config); + } + break; + default: + SU2_MPI::Error("Invalid upwind scheme or not implemented.\n Currently, only UDS is available for pressure based incompressible flows.", CURRENT_FUNCTION); + break; + } } } break; @@ -2095,6 +2127,21 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } } + /*--- Solver definition for the poisson/pressure correction problem ---*/ + if (poisson) { + /*--- Pressure correction (Poisson) equation ---*/ + numerics[MESH_0][POISSON_SOL][visc_term] = new CAvgGrad_Heat(nDim, config, true); + + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][POISSON_SOL][visc_term] = new CAvgGrad_Heat(nDim, config, false); + + /*--- Assign the convective boundary term as well to account for flow BCs as well --*/ + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { + numerics[iMGlevel][POISSON_SOL][visc_bound_term] = new CAvgGrad_Heat(nDim, config, false); + + } + } + /*--- Solver definition for the radiation model problem ---*/ if (config->AddRadiation()) { diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index 468e5cbb231..66a9927a43c 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -87,7 +87,39 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe const bool fmg_cfl_ramp = integration[val_iZone][val_iInst][FLOW_SOL]->GetFullMG_CFLRamp(); /*--- If the flow integration is not fully coupled, run the various single grid integrations. ---*/ + CommonAuxiliarySolvers(output, integration, geometry, solver, numerics, config, surface_movement, + grid_movement, FFDBox, val_iZone, val_iInst, main_solver, frozen_visc); + /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. + The Full-MG startup owns the CFL while it ramps, so leave it alone until then. ---*/ + SU2_OMP_PARALLEL + if (!disc_adj && config[val_iZone]->GetFinestMesh() == MESH_0 && !fmg_cfl_ramp) { + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], + solver[val_iZone][val_iInst], config[val_iZone]); + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->IdentifySolutionOutliers(config[val_iZone], InnerIter); + } + END_SU2_OMP_PARALLEL + + /*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/ + + if ((config[val_iZone]->GetGrid_Movement()) && (config[val_iZone]->GetAeroelastic_Simulation()) && unsteady) { + SetGrid_Movement(geometry[val_iZone][val_iInst], surface_movement[val_iZone], grid_movement[val_iZone][val_iInst], + solver[val_iZone][val_iInst], config[val_iZone], InnerIter, TimeIter); + + /*--- Apply a Wind Gust ---*/ + + if (config[val_iZone]->GetWind_Gust()) { + if (InnerIter % config[val_iZone]->GetAeroelasticIter() == 0 && InnerIter != 0) + SetWind_GustField(config[val_iZone], geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst]); + } + } +} + +void CFluidIteration::CommonAuxiliarySolvers(COutput* output, CIntegration**** integration, CGeometry**** geometry, + CSolver***** solver, CNumerics****** numerics, CConfig** config, + CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, MAIN_SOLVER main_solver, bool frozen_visc) { + if (config[val_iZone]->GetKind_Turb_Model() != TURB_MODEL::NONE && !frozen_visc) { /*--- Solve transition model ---*/ @@ -133,30 +165,7 @@ void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGe integration[val_iZone][val_iInst][RAD_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_RADIATION_SYS, val_iZone, val_iInst); } - - /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. - The Full-MG startup owns the CFL while it ramps, so leave it alone until then. ---*/ - SU2_OMP_PARALLEL - if (!disc_adj && config[val_iZone]->GetFinestMesh() == MESH_0 && !fmg_cfl_ramp) { - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], - solver[val_iZone][val_iInst], config[val_iZone]); - solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->IdentifySolutionOutliers(config[val_iZone], InnerIter); - } - END_SU2_OMP_PARALLEL - - /*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/ - - if ((config[val_iZone]->GetGrid_Movement()) && (config[val_iZone]->GetAeroelastic_Simulation()) && unsteady) { - SetGrid_Movement(geometry[val_iZone][val_iInst], surface_movement[val_iZone], grid_movement[val_iZone][val_iInst], - solver[val_iZone][val_iInst], config[val_iZone], InnerIter, TimeIter); - - /*--- Apply a Wind Gust ---*/ - - if (config[val_iZone]->GetWind_Gust()) { - if (InnerIter % config[val_iZone]->GetAeroelasticIter() == 0 && InnerIter != 0) - SetWind_GustField(config[val_iZone], geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst]); - } - } + } void CFluidIteration::Update(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, diff --git a/SU2_CFD/src/iteration/CIterationFactory.cpp b/SU2_CFD/src/iteration/CIterationFactory.cpp index 761ceca98f8..6d738b62b28 100644 --- a/SU2_CFD/src/iteration/CIterationFactory.cpp +++ b/SU2_CFD/src/iteration/CIterationFactory.cpp @@ -32,6 +32,7 @@ #include "../../include/iteration/CDiscAdjFluidIteration.hpp" #include "../../include/iteration/CDiscAdjHeatIteration.hpp" #include "../../include/iteration/CFluidIteration.hpp" +#include "../../include/iteration/CPBFluidIteration.hpp" #include "../../include/iteration/CFEMFluidIteration.hpp" #include "../../include/iteration/CTurboIteration.hpp" #include "../../include/iteration/CHeatIteration.hpp" @@ -57,6 +58,11 @@ CIteration* CIterationFactory::CreateIteration(MAIN_SOLVER kindSolver, const CCo iteration = new CTurboIteration(config); } + else if (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED) { + if (rank == MASTER_NODE) + cout << "Pressure based Euler/Navier-Stokes/RANS fluid iteration." << endl; + iteration = new CPBFluidIteration(config); + } else{ if (rank == MASTER_NODE) cout << "Euler/Navier-Stokes/RANS fluid iteration." << endl; @@ -113,7 +119,7 @@ CIteration* CIterationFactory::CreateIteration(MAIN_SOLVER kindSolver, const CCo iteration = new CDiscAdjHeatIteration(config); break; - case MAIN_SOLVER::NONE: case MAIN_SOLVER::TEMPLATE_SOLVER: case MAIN_SOLVER::MULTIPHYSICS: + case MAIN_SOLVER::NONE: case MAIN_SOLVER::TEMPLATE_SOLVER: case MAIN_SOLVER::MULTIPHYSICS: case MAIN_SOLVER::POISSON_EQUATION: SU2_MPI::Error("No iteration found for specified solver.", CURRENT_FUNCTION); break; } diff --git a/SU2_CFD/src/iteration/CPBFluidIteration.cpp b/SU2_CFD/src/iteration/CPBFluidIteration.cpp new file mode 100644 index 00000000000..bd4f886ec2f --- /dev/null +++ b/SU2_CFD/src/iteration/CPBFluidIteration.cpp @@ -0,0 +1,102 @@ +/*! + * \file CPBFluidIteration.cpp + * \brief Main subroutines used by SU2_CFD + * \author F. Palacios, T. Economon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/iteration/CPBFluidIteration.hpp" +#include "../../include/output/COutput.hpp" +#include "../../include/integration/CIntegration.hpp" + +void CPBFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, + CSolver***** solver, CNumerics****** numerics, CConfig** config, + CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + SU2_ZONE_SCOPED + + const bool frozen_visc = (config[val_iZone]->GetContinuous_Adjoint() && config[val_iZone]->GetFrozen_Visc_Cont()) || + (config[val_iZone]->GetDiscrete_Adjoint() && config[val_iZone]->GetFrozen_Visc_Disc()); + const bool disc_adj = (config[val_iZone]->GetDiscrete_Adjoint()); + const bool periodic = (config[val_iZone]->GetnMarker_Periodic() > 0); + + const unsigned short nCorrections = config[val_iZone]->GetSIMPLE_Options().nCorrections_PISO; + + /*--- Solve the Euler, Navier-Stokes, RANS equations. ---*/ + + const auto main_solver = config[val_iZone]->GetKind_Solver(); + config[val_iZone]->SetGlobalParam(main_solver, RUNTIME_FLOW_SYS); + + /*--- Solve the momentum equations (to find the predicted velocity u*). ---*/ + + integration[val_iZone][val_iInst][FLOW_SOL]->MultiGrid_Iteration(geometry, solver, numerics, config, RUNTIME_FLOW_SYS, + val_iZone, val_iInst); + + /*--- The momentum coefficients (resulting from the flow solution) are set only once + at the start of the corrections. These coefficients make up the entirety of the coefficient + matrix (Jacobian) used by the Poisson solver. Currently the matrix is redefined each correction + but as the coefficients are frozen this doesnt/shouldnt change the matrix at all. ---*/ + + solver[val_iZone][val_iInst][MESH_0][POISSON_SOL]->SetMomCoeff(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone], periodic, MESH_0); + + /*--- Compute the mass fluxes at the cell edges based on Rhie-Chow interpolation ---*/ + + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->ComputeEdgeMassFluxesRhieChow(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); + + /*--- Solve the pressure poisson (correction) equation ---*/ + + config[val_iZone]->SetGlobalParam(MAIN_SOLVER::POISSON_EQUATION, RUNTIME_POISSON_SYS); + + for (unsigned short iCorrection = 0; iCorrection < nCorrections; ++iCorrection) { + + /*--- For later corrections (PISO) the pressure equation has an additional div(H(u')/A_p) term on the right side, u' is computed in the last correction routine. ---*/ + + if (iCorrection > 0) solver[val_iZone][val_iInst][MESH_0][POISSON_SOL]->ComputeHbyA(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone], MESH_0); + + /*--- Solve the pressure Poisson equation to find p' i.e. div(V/A_p * grad(p')) = div rhou*} ---*/ + + integration[val_iZone][val_iInst][POISSON_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config, RUNTIME_POISSON_SYS, + val_iZone, val_iInst); + + /*--- The velocity and pressure are corrected based on the solution to the Poisson problem i.e. p* = p + p' and rhou** = rhou* - V/Ap * p' ---*/ + + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->ApplyPressureVelocityCorrection(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); + + } + + /*--- Pressure-based algorithm finished, now run auxiliary solvers ---*/ + + /*--- If the flow integration is not fully coupled, run the various single grid integrations. ---*/ + CommonAuxiliarySolvers(output, integration, geometry, solver, numerics, config, surface_movement, + grid_movement, FFDBox, val_iZone, val_iInst, main_solver, frozen_visc); + + /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/ + + if ((config[val_iZone]->GetCFL_Adapt() == YES) && (!disc_adj)) { + SU2_OMP_PARALLEL + solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], + solver[val_iZone][val_iInst], config[val_iZone]); + END_SU2_OMP_PARALLEL + } + +} diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index d4db53843e9..4462b6b4111 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -82,6 +82,7 @@ su2_cfd_src += files(['variables/CIncNSVariable.cpp', 'variables/CAdjTurbVariable.cpp', 'variables/CFlowVariable.cpp', 'variables/CIncEulerVariable.cpp', + 'variables/CPoissonVariable.cpp', 'variables/CEulerVariable.cpp', 'variables/CNEMOEulerVariable.cpp', 'variables/CNEMONSVariable.cpp', @@ -105,6 +106,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CHeatSolver.cpp', 'solvers/CIncEulerSolver.cpp', 'solvers/CIncNSSolver.cpp', + 'solvers/CPoissonSolver.cpp', 'solvers/CMeshSolver.cpp', 'solvers/CNEMOEulerSolver.cpp', 'solvers/CNEMONSSolver.cpp', @@ -129,6 +131,7 @@ su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/flow/convection/hllc.cpp', 'numerics/flow/convection/ausm_slau.cpp', 'numerics/flow/convection/centered.cpp', + 'numerics/flow/convection/pressure_based.cpp', 'numerics/flow/flow_diffusion.cpp', 'numerics/flow/flow_sources.cpp', 'numerics/NEMO/CNEMONumerics.cpp', @@ -185,6 +188,7 @@ su2_cfd_src += files(['iteration/CIteration.cpp', 'iteration/CFEAIteration.cpp', 'iteration/CFEMFluidIteration.cpp', 'iteration/CFluidIteration.cpp', + 'iteration/CPBFluidIteration.cpp', 'iteration/CHeatIteration.cpp', 'iteration/CTurboIteration.cpp']) diff --git a/SU2_CFD/src/numerics/flow/convection/pressure_based.cpp b/SU2_CFD/src/numerics/flow/convection/pressure_based.cpp new file mode 100644 index 00000000000..f6fe14715af --- /dev/null +++ b/SU2_CFD/src/numerics/flow/convection/pressure_based.cpp @@ -0,0 +1,198 @@ +/*! + * \file pressure_based.cpp + * \brief Implementations of fluxes for pressure-based solvers. + * \author T. Aalbers + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../../../include/numerics/flow/convection/pressure_based.hpp" + +CPBConvection_Base::CPBConvection_Base(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + dynamic_grid = config->GetDynamic_Grid(); + energy = config->GetEnergy_Equation(); + variable_density = (config->GetVariable_Density_Model()); + + AdvectedVelocity = new su2double [MAXNDIM]; + Flux = new su2double [nVar]; + Jacobian_i = new su2double* [nVar]; + Jacobian_j = new su2double* [nVar]; + + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nVar]; + Jacobian_j[iVar] = new su2double [nVar]; + } +} + +CPBConvection_Base::~CPBConvection_Base(void) { + + delete [] AdvectedVelocity; + delete [] Flux; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + } + + delete [] Jacobian_i; + delete [] Jacobian_j; + +} + +CNumerics::ResidualType<> CPBConvection_Base::ComputeResidual(const CConfig *config) { + + /*--- Primitive variables at point i and j ---*/ + + Pressure_i = V_i[0]; Pressure_j = V_j[0]; + DensityInc_i = V_i[nDim+2]; DensityInc_j = V_j[nDim+2]; + Enthalpy_i = V_i[nDim+3]; Enthalpy_j = V_j[nDim+3]; + MeanPressure = 0.5 * (Pressure_i + Pressure_j); + MeanDensity = 0.5 * (DensityInc_i + DensityInc_j); + + /*--- Find the velocity that is advected ---*/ + + ComputeAdvectedQuantities(); + + /*--- Set the flux vector. ---*/ + + Flux[0] = MassFlux; + for (iDim = 0; iDim < nDim; ++iDim) + Flux[1+iDim] = MassFlux * AdvectedVelocity[iDim]; + Flux[nDim+1] = MassFlux * AdvectedEnthalpy; + + /*--- Find Jacobian ---*/ + + if (implicit) { + + for (jVar = 0; jVar < nVar; jVar++) + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_j[iVar][jVar] = 0.0; + } + + /*--- We need the derivative of the equation of state to build the + preconditioning matrix. For now, the only option is the ideal gas + law, but in the future, dRhodT should be in the fluid model. ---*/ + + dRhodh_i = 0.0; dRhodh_j = 0.0; + if (variable_density) { + Temperature_i = V_i[nDim+1]; Temperature_j = V_j[nDim+1]; + Cp_i = V_i[nDim+8]; Cp_j = V_j[nDim+8]; + + dRhodh_i = -DensityInc_i / (Temperature_i * Cp_i); + dRhodh_j = -DensityInc_j / (Temperature_j * Cp_j); + } + + ComputeJacobianWeights(); + ComputeJacobian(DensityInc_i, &V_i[1], Enthalpy_i, dRhodh_i, weight_jacobian_i, Jacobian_i); + ComputeJacobian(DensityInc_j, &V_j[1], Enthalpy_j, dRhodh_j, weight_jacobian_j, Jacobian_j); + + } + + /*--- Remove energy contributions if we aren't solving the energy equation. ---*/ + + if (!energy) { + Flux[nDim+1] = 0.0; + if (implicit) { + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar][nDim+1] = 0.0; + Jacobian_j[iVar][nDim+1] = 0.0; + + Jacobian_i[nDim+1][iVar] = 0.0; + Jacobian_j[nDim+1][iVar] = 0.0; + } + } + } + + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); +} + + +void CPBConvection_Base::ComputeJacobian(su2double val_density, const su2double *val_velocity, + su2double val_enthalpy, su2double val_dRhodh, + su2double val_scale, su2double **val_Proj_Jac_Tensor) { + + su2double proj_vel = MassFlux / MeanDensity; + + /*--- Fill continuity parts ---*/ + + val_Proj_Jac_Tensor[0][0] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + val_Proj_Jac_Tensor[0][iDim+1] = val_scale*(Normal[iDim] * (val_density)); + } + + /*--- Fill momentum parts ---*/ + + for (jDim = 0; jDim < nDim; jDim++) { + for (iDim = 0; iDim < nDim; iDim++) { + val_Proj_Jac_Tensor[iDim+1][jDim+1] = val_scale * val_density * (val_velocity[iDim] * Normal[jDim] + proj_vel * delta[iDim][jDim]); + } + } + + /*--- Fill enthalpy parts ---*/ + + val_Proj_Jac_Tensor[nDim+1][0] = 0.0; + val_Proj_Jac_Tensor[0][nDim+1] = val_scale * ((val_dRhodh) * proj_vel); + val_Proj_Jac_Tensor[nDim+1][nDim+1] = val_scale*(((val_enthalpy)*(val_dRhodh) + (val_density))*proj_vel); + for (iDim = 0; iDim < nDim; iDim++) { + val_Proj_Jac_Tensor[nDim+1][iDim+1] = val_scale * ((val_enthalpy) * Normal[iDim] * (val_density)); + val_Proj_Jac_Tensor[iDim+1][nDim+1] = val_scale*((val_dRhodh) * val_velocity[iDim] * proj_vel); + } + +} + +void CPBConvection_Central::ComputeAdvectedQuantities() { + + for (iDim = 0; iDim < nDim; iDim++) + AdvectedVelocity[iDim] = 0.5 * (V_i[iDim+1] + V_j[iDim+1]); + + AdvectedEnthalpy = 0.5 * (Enthalpy_i + Enthalpy_j); + +} + +void CPBConvection_Central::ComputeJacobianWeights() { + + weight_jacobian_i = weight_jacobian_j = 0.5; + +} + +void CPBConvection_Upwind::ComputeAdvectedQuantities() { + + bool Upw_i = (MassFlux>0); + + for (iDim = 0; iDim < nDim; iDim++) + AdvectedVelocity[iDim] = Upw_i ? V_i[iDim+1] : V_j[iDim+1]; + + AdvectedEnthalpy = (Upw_i) ? Enthalpy_i : Enthalpy_j; + +} + +void CPBConvection_Upwind::ComputeJacobianWeights() { + + bool Upw_i = (MassFlux>0); + + weight_jacobian_i = static_cast(Upw_i); + weight_jacobian_j = static_cast(!Upw_i); + +} diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 23e5f066f15..bd25cc9712b 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -40,6 +40,8 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); flamelet = (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET); + pressure_based = (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED); + streamwisePeriodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE); streamwisePeriodic_temperature = config->GetStreamwise_Periodic_Temperature(); @@ -91,8 +93,12 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp restartFilename = config->GetRestart_FileName(); /*--- Set the default convergence field --- */ + /*--- The default field for the pressure-based solver is velocity, as the pressure is only + solved as a correction and it is thus possible to run the pb solver without corrections. to + ensure the default residual is always defined we use velocity instead ---*/ - if (convFields.empty()) convFields.emplace_back("RMS_PRESSURE"); + if (convFields.empty() && !pressure_based) convFields.emplace_back("RMS_PRESSURE"); + if (convFields.empty() && pressure_based) convFields.emplace_back("RMS_VELOCITY-X"); } @@ -111,7 +117,7 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ if (weakly_coupled_heat) AddHistoryOutput("RMS_TEMPERATURE", "rms[T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the temperature.", HistoryFieldType::RESIDUAL); /// DESCRIPTION: Root-mean square residual of the enthalpy. if (heat) AddHistoryOutput("RMS_ENTHALPY", "rms[h]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the enthalpy.", HistoryFieldType::RESIDUAL); - + AddHistoryOutputFields_ScalarRMS_RES(config); /// DESCRIPTION: Root-mean square residual of the radiative energy (P1 model). @@ -168,6 +174,10 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){ /// DESCRIPTION: Linear solver iterations AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); + if (pressure_based) { + AddHistoryOutput("LINSOL_POISSON_ITER", "PoissonSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the poisson solver."); + AddHistoryOutput("LINSOL_POISSON_RESIDUAL", "PoissonSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the poisson solver."); + } AddHistoryOutputFieldsScalarLinsol(config); AddHistoryOutput("MIN_DELTA_TIME", "Min DT", ScreenOutputFormat::SCIENTIFIC, "CFL_NUMBER", "Current minimum local time step"); @@ -209,8 +219,13 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv CSolver* heat_solver = solver[HEAT_SOL]; CSolver* rad_solver = solver[RAD_SOL]; CSolver* mesh_solver = solver[MESH_SOL]; + CSolver* poisson_solver = solver[POISSON_SOL]; - SetHistoryOutputValue("RMS_PRESSURE", log10(flow_solver->GetRes_RMS(0))); + if (pressure_based) { + SetHistoryOutputValue("RMS_PRESSURE", log10(poisson_solver->GetRes_RMS(0))); + } else { + SetHistoryOutputValue("RMS_PRESSURE", log10(flow_solver->GetRes_RMS(0))); + } SetHistoryOutputValue("RMS_VELOCITY-X", log10(flow_solver->GetRes_RMS(1))); SetHistoryOutputValue("RMS_VELOCITY-Y", log10(flow_solver->GetRes_RMS(2))); if (nDim == 3) SetHistoryOutputValue("RMS_VELOCITY-Z", log10(flow_solver->GetRes_RMS(3))); @@ -218,7 +233,11 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv if (config->AddRadiation()) SetHistoryOutputValue("RMS_RAD_ENERGY", log10(rad_solver->GetRes_RMS(0))); - SetHistoryOutputValue("MAX_PRESSURE", log10(flow_solver->GetRes_Max(0))); + if (pressure_based) { + SetHistoryOutputValue("MAX_PRESSURE", log10(poisson_solver->GetRes_Max(0))); + } else { + SetHistoryOutputValue("MAX_PRESSURE", log10(flow_solver->GetRes_Max(0))); + } SetHistoryOutputValue("MAX_VELOCITY-X", log10(flow_solver->GetRes_Max(1))); SetHistoryOutputValue("MAX_VELOCITY-Y", log10(flow_solver->GetRes_Max(2))); if (nDim == 3) SetHistoryOutputValue("MAX_VELOCITY-Z", log10(flow_solver->GetRes_Max(3))); @@ -251,6 +270,10 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv SetHistoryOutputValue("LINSOL_ITER", flow_solver->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(flow_solver->GetResLinSolver())); + if (pressure_based) { + SetHistoryOutputValue("LINSOL_POISSON_ITER",poisson_solver->GetIterLinSolver()); + SetHistoryOutputValue("LINSOL_POISSON_RESIDUAL",log10(poisson_solver->GetResLinSolver())); + } if (config->GetDeform_Mesh()){ SetHistoryOutputValue("DEFORM_MIN_VOLUME", mesh_solver->GetMinimum_Volume()); @@ -447,7 +470,11 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint)); } - SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0)); + if (pressure_based) { + SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[POISSON_SOL]->LinSysRes(iPoint, 0)); + } else { + SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0)); + } SetVolumeOutputValue("RES_VELOCITY-X", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 1)); SetVolumeOutputValue("RES_VELOCITY-Y", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 2)); if (nDim == 3) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index acae4cbde05..32c265e5565 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -165,7 +165,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, false, true); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, true); } else { if (rank == MASTER_NODE) diff --git a/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp b/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp index fb26eee6661..a439f913c83 100644 --- a/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp +++ b/SU2_CFD/src/solvers/CGradientSmoothingSolver.cpp @@ -107,18 +107,20 @@ CGradientSmoothingSolver::CGradientSmoothingSolver(CGeometry *geometry, CConfig } /*--- initializations for linear equation systems ---*/ + std::optional override_prec = + config->GetSmoothGradient() ? std::optional{config->GetKind_Grad_Linear_Solver_Prec()} : std::nullopt; if ( !config->GetSmoothOnSurface() ) { nVar = config->GetSmoothSepDim() ? 1 : nDim; LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, false, geometry, config, false, true); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, false, geometry, config, false, false, override_prec); } else { if (config->GetSobMode() == ENUM_SOBOLEV_MODUS::PARAM_LEVEL_COMPLETE) { - Jacobian.Initialize(nPoint, nPointDomain, nDim, nDim, false, geometry, config, false , true); + Jacobian.Initialize(nPoint, nPointDomain, nDim, nDim, false, geometry, config, false, false, override_prec); } else { LinSysSol.Initialize(nPoint, nPointDomain, 1, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, 1, 0.0); - Jacobian.Initialize(nPoint, nPointDomain, 1, 1, false, geometry, config, false, true); + Jacobian.Initialize(nPoint, nPointDomain, 1, 1, false, geometry, config, false, false, override_prec); } visited.resize(geometry->GetnPoint(), false); } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 41580c3062a..e37a062928e 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -39,7 +39,8 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, const bool navier_stokes) : - CFVMFlowSolverBase(*geometry, *config) { + CFVMFlowSolverBase(*geometry, *config), + pressure_based(config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED) { SU2_ZONE_SCOPED /*--- Based on the navier_stokes boolean, determine if this constructor is @@ -155,7 +156,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, false, true); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, true); } else { if (rank == MASTER_NODE) @@ -216,8 +217,30 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned CommunicateInitialState(geometry, config); /*--- Sizing edge mass flux array ---*/ - if (config->GetBounded_Scalar()) + if (config->GetBounded_Scalar() || pressure_based) EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0); + + /*--- Pressure based solver specific allocations ---*/ + if (pressure_based) { + + /*--- Initialize the edge mass flux array ---*/ + + for (unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + + EdgeMassFluxes[iEdge] = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + EdgeMassFluxes[iEdge] += Density_Inf * Velocity_Inf[iDim] * geometry->edges->GetNormal(iEdge)[iDim]; + + } + + /*--- Allocate corrections and relaxation ---*/ + + pressureCorrection.resize(nPointDomain) = su2double(0.0); + momentumCorrection.resize(nPointDomain,nDim) = su2double(0.0); + EdgeMassFluxCorrection.resize(geometry->GetnEdge()) = su2double(0.0); + alpha_p.resize(nPointDomain) = su2double(1.0); + + } /*--- Add the solver name. ---*/ SolverName = "INC.FLOW"; @@ -1003,7 +1026,7 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ /*--- Update the beta value based on the maximum velocity. ---*/ - SetBeta_Parameter(geometry, solver_container, config, iMesh); + if (!pressure_based) SetBeta_Parameter(geometry, solver_container, config, iMesh); /*--- Update the pressure range in the domain for target outflow mass flow rate. ---*/ @@ -1015,6 +1038,14 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ SU2_OMP_SAFE_GLOBAL_ACCESS(GetOutlet_Properties(geometry, config, iMesh, Output);) } + /*--- Reset flag for strong BCs. ---*/ + if (pressure_based) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + nodes->ResetStrongBC(iPoint); + END_SU2_OMP_FOR + } + /*--- Initialize the Jacobian matrix and residual, not needed for the reducer strategy * as we set blocks (including diagonal ones) and completely overwrite. ---*/ @@ -1124,18 +1155,6 @@ void CIncEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_contain unsigned short iMesh, unsigned long Iteration) { SU2_ZONE_SCOPED - /*--- Define an object to compute the speed of sound. ---*/ - struct SoundSpeed { - FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { - return sqrt(0.5 * (nodes.GetBetaInc2(iPoint) + nodes.GetBetaInc2(jPoint))); - } - - FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { - return sqrt(nodes.GetBetaInc2(iPoint)); - } - - } soundSpeed; - /*--- Define an object to compute the viscous eigenvalue. ---*/ struct LambdaVisc { const bool energy; @@ -1169,9 +1188,38 @@ void CIncEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_contain } lambdaVisc(config->GetEnergy_Equation()); - /*--- Now instantiate the generic implementation with the two functors above. ---*/ + if (pressure_based) { + /* Define an object to compute the speed of sound, as the speed of sound is theoretically infinite, + this makes no sense. However to be able to reuse the time step routine we artificially define the speed of sound + to be zero such that a regular advective time step is computed */ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint = 0) const { + return 0.0; + } + + } soundSpeed; - SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); + /*--- Now instantiate the generic implementation with the two functors above. ---*/ + + SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); + + } else { + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return sqrt(0.5 * (nodes.GetBetaInc2(iPoint) + nodes.GetBetaInc2(jPoint))); + } + + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { + return sqrt(nodes.GetBetaInc2(iPoint)); + } + + } soundSpeed; + + /*--- Now instantiate the generic implementation with the two functors above. ---*/ + + SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); + } } @@ -1240,6 +1288,10 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(jPoint)); } + /*--- Set the edge mass flux ---*/ + + if (pressure_based) numerics->SetMassFlux(EdgeMassFluxes[iEdge]); + /*--- Compute residuals, and Jacobians ---*/ auto conv_residual = numerics->ComputeResidual(config); @@ -1413,6 +1465,10 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont } + /*--- Set the edge mass flux ---*/ + + if (pressure_based) numerics->SetMassFlux(EdgeMassFluxes[iEdge]); + /*--- Compute the residual ---*/ auto conv_residual = numerics->ComputeResidual(config); @@ -1477,6 +1533,20 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont AD::StartNoSharedReading(); + if (pressure_based) { + + /*--- Add pressure source term (V * gradp) ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + for (unsigned short iDim = 0; iDim < nDim; iDim++) + LinSysRes(iPoint, iDim + 1) += geometry->nodes->GetVolume(iPoint) * nodes->GetGradient_Primitive(iPoint,prim_idx.Pressure(),iDim); + + } + END_SU2_OMP_FOR + } + if (body_force) { /*--- Loop over all points ---*/ @@ -1984,6 +2054,12 @@ template FORCEINLINE void CIncEulerSolver::Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep) { SU2_ZONE_SCOPED + + if (pressure_based) { + CFVMFlowSolverBase::Explicit_Iteration(geometry, solver_container, config, iRKStep); + return; + } + struct Precond { const CIncEulerSolver* solver; su2activematrix matrix; @@ -2033,10 +2109,12 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C struct IncPrec { const CIncEulerSolver* solver; - const bool active = true; + const bool active; su2activematrix matrix; - IncPrec(const CIncEulerSolver* s, unsigned short nVar) : solver(s) { matrix.resize(nVar,nVar); } + IncPrec(const CIncEulerSolver* s, unsigned short nVar) : solver(s), active(!s->pressure_based) { + matrix.resize(nVar,nVar); + } FORCEINLINE const su2activematrix& operator() (const CConfig* config, unsigned long iPoint, su2double delta) { solver->SetPreconditioner(config, iPoint, delta, matrix); @@ -2046,6 +2124,17 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C } precond(this, nVar); PrepareImplicitIteration_impl(precond, geometry, config); + + /*--- Delete pressure rows for segregated solver type. ---*/ + if (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED) { + SU2_OMP_FOR_(schedule(static,omp_chunk_size) SU2_NOWAIT) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + Jacobian.DeleteValsRowi(iPoint, 0); + LinSysRes(iPoint,0) = 0.0; + LinSysSol(iPoint,0) = 0.0; + } + END_SU2_OMP_FOR + } } void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_container, @@ -2292,7 +2381,6 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain V_infty[prim_idx.Pressure()] = GetPressure_Inf(); /*--- Dirichlet condition for temperature at far-field (if energy is active). ---*/ - V_infty[prim_idx.Temperature()] = GetTemperature_Inf(); /*-- Enthalpy at far-field. ---*/ @@ -2308,7 +2396,7 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain /*--- Beta coefficient stored at the node ---*/ - V_infty[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); + if (!pressure_based) V_infty[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); /*--- Cp is needed for Temperature equation. ---*/ @@ -2322,22 +2410,75 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Compute the convective residual using an upwind scheme ---*/ + if (pressure_based) { - auto residual = conv_numerics->ComputeResidual(config); + /*--- Decide if the boundary should be an inlet or an outlet ---*/ - /*--- Update residual value ---*/ + su2double Face_Flux = 0.0; + if (dynamic_grid) + for (iDim = 0; iDim < nDim; iDim++) + Face_Flux += nodes->GetDensity(iPoint)*(V_domain[iDim+1]-geometry->nodes->GetGridVel(iPoint)[iDim])*Normal[iDim]; + else + for (iDim = 0; iDim < nDim; iDim++) + Face_Flux += nodes->GetDensity(iPoint)*V_domain[iDim+1]*Normal[iDim]; - LinSysRes.AddBlock(iPoint, residual); + bool inflow = false; + if ((Face_Flux < 0.0) && (fabs(Face_Flux) > EPS)) inflow = true; - /*--- Convective Jacobian contribution for implicit integration ---*/ + if (inflow) { - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + /*--- Set this face as an inlet via a strong BC. ---*/ + + LinSysRes.SetBlock_Zero(iPoint); + + /*--- Mark as a strong BC which is important for deciding if a velocity correction should be applied ---*/ + + nodes->SetStrongBC(iPoint); + + if (implicit) + for (iDim = 0; iDim < nDim; iDim++) + Jacobian.DeleteValsRowi(iPoint, iDim+1); + + } else { + + /*--- Set the edge mass flux ---*/ + + conv_numerics->SetMassFlux(Face_Flux); + + /*--- Compute the residual using an upwind scheme ---*/ + + conv_numerics->SetPrimitive(V_domain, V_domain); + + auto residual = conv_numerics->ComputeResidual(config); + + LinSysRes.AddBlock(iPoint, residual); + nodes->SetSolution(iPoint, 0, GetPressure_Inf()); + nodes->SetPressure(iPoint); + + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + + } + } else { + + /*--- Compute the convective residual using an upwind scheme ---*/ + + auto residual = conv_numerics->ComputeResidual(config); + + /*--- Update residual value ---*/ + + LinSysRes.AddBlock(iPoint, residual); + + /*--- Convective Jacobian contribution for implicit integration ---*/ + + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + + } /*--- Viscous residual contribution ---*/ - if (!viscous || energy_multicomponent) continue; + if (!viscous || energy_multicomponent || pressure_based) continue; /*--- Set transport properties at infinity. ---*/ @@ -2481,6 +2622,9 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, case INLET_TYPE::PRESSURE_INLET: + if (pressure_based) + SU2_MPI::Error("Pressure Inlet is currently an unsupported INC_INLET_TYPE for pressure based solver.", CURRENT_FUNCTION); + /*--- Retrieve the specified total pressure for the inlet. ---*/ P_total = Inlet_Ptotal[val_marker][iVertex]/config->GetPressure_Ref(); @@ -2558,96 +2702,111 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[prim_idx.Pressure()] = nodes->GetPressure(iPoint); } - /*-- Enthalpy is needed for energy equation. ---*/ - const su2double* scalar_inlet = nullptr; - if (species_model) scalar_inlet = config->GetInlet_SpeciesVal(config->GetMarker_All_TagBound(val_marker)); - CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel(); - auxFluidModel->SetTDState_T(V_inlet[prim_idx.Temperature()], scalar_inlet); + if (pressure_based) { - /*--- For the flamelet model with FLOW_MARKERS enthalpy BC, we obtain the inlet enthalpy - from the flamelet species solver With SPECIES_MARKERS, the enthalpy in MARKER_INLET_SPECIES - is used directly. ---*/ - if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET && - config->GetFlamelet_Enthalpy_BC() == FLAMELET_ENTHALPY_BC::FLOW_MARKERS) - V_inlet[prim_idx.Enthalpy()] = nodes->GetEnthalpy(iPoint); - else - V_inlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy(); + /*--- Directly overwrite the velocity at the boundary nodes as a dirichlet boundary condition ---*/ - /*--- Access density at the node. This is either constant by - construction, or will be set fixed implicitly by the temperature - and equation of state. ---*/ + nodes->SetVelocity_Old(iPoint,V_inlet+1); - V_inlet[prim_idx.Density()] = nodes->GetDensity(iPoint); + LinSysRes.SetBlock_Zero(iPoint); - /*--- Beta coefficient from the config file ---*/ + if (implicit) + for (iDim = 0; iDim < nDim; iDim++) + Jacobian.DeleteValsRowi(iPoint, iDim+1); - V_inlet[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); + } else { - /*--- Cp is needed for Temperature equation. ---*/ + /*-- Enthalpy is needed for energy equation. ---*/ + const su2double* scalar_inlet = nullptr; + if (species_model) scalar_inlet = config->GetInlet_SpeciesVal(config->GetMarker_All_TagBound(val_marker)); + CFluidModel* auxFluidModel = solver_container[FLOW_SOL]->GetFluidModel(); + auxFluidModel->SetTDState_T(V_inlet[prim_idx.Temperature()], scalar_inlet); + + /*--- For the flamelet model with FLOW_MARKERS enthalpy BC, we obtain the inlet enthalpy + from the flamelet species solver With SPECIES_MARKERS, the enthalpy in MARKER_INLET_SPECIES + is used directly. ---*/ + if (config->GetKind_Species_Model() == SPECIES_MODEL::FLAMELET && + config->GetFlamelet_Enthalpy_BC() == FLAMELET_ENTHALPY_BC::FLOW_MARKERS) + V_inlet[prim_idx.Enthalpy()] = nodes->GetEnthalpy(iPoint); + else + V_inlet[prim_idx.Enthalpy()] = auxFluidModel->GetEnthalpy(); + + /*--- Access density at the node. This is either constant by + construction, or will be set fixed implicitly by the temperature + and equation of state. ---*/ - V_inlet[prim_idx.CpTotal()] = nodes->GetSpecificHeatCp(iPoint); + V_inlet[prim_idx.Density()] = nodes->GetDensity(iPoint); - /*--- Set various quantities in the solver class ---*/ + /*--- Beta coefficient from the config file ---*/ - conv_numerics->SetPrimitive(V_domain, V_inlet); + V_inlet[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + /*--- Cp is needed for Temperature equation. ---*/ - /*--- Compute the residual using an upwind scheme ---*/ + V_inlet[prim_idx.CpTotal()] = nodes->GetSpecificHeatCp(iPoint); - auto residual = conv_numerics->ComputeResidual(config); + /*--- Set various quantities in the solver class ---*/ - /*--- Update residual value ---*/ + conv_numerics->SetPrimitive(V_domain, V_inlet); - LinSysRes.AddBlock(iPoint, residual); + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Compute the residual using an upwind scheme ---*/ - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + auto residual = conv_numerics->ComputeResidual(config); - /*--- Viscous contribution, commented out because serious convergence problems ---*/ + /*--- Update residual value ---*/ - if (!viscous || energy_multicomponent) continue; + LinSysRes.AddBlock(iPoint, residual); - /*--- Set transport properties at the inlet ---*/ + /*--- Jacobian contribution for implicit integration ---*/ - V_inlet[prim_idx.LaminarViscosity()] = nodes->GetLaminarViscosity(iPoint); - V_inlet[prim_idx.EddyViscosity()] = nodes->GetEddyViscosity(iPoint); - V_inlet[prim_idx.ThermalConductivity()] = nodes->GetThermalConductivity(iPoint); + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - /*--- Set the normal vector and the coordinates ---*/ + /*--- Viscous contribution, commented out because serious convergence problems ---*/ - visc_numerics->SetNormal(Normal); - su2double Coord_Reflected[MAXNDIM]; - GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), - geometry->nodes->GetCoord(iPoint), Coord_Reflected); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected); + if (!viscous || energy_multicomponent) continue; - /*--- Primitive variables, and gradient ---*/ + /*--- Set transport properties at the inlet ---*/ - visc_numerics->SetPrimitive(V_domain, V_inlet); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - nodes->GetGradient_Primitive(iPoint)); + V_inlet[prim_idx.LaminarViscosity()] = nodes->GetLaminarViscosity(iPoint); + V_inlet[prim_idx.EddyViscosity()] = nodes->GetEddyViscosity(iPoint); + V_inlet[prim_idx.ThermalConductivity()] = nodes->GetThermalConductivity(iPoint); - /*--- Turbulent kinetic energy ---*/ + /*--- Set the normal vector and the coordinates ---*/ - if (config->GetKind_Turb_Model() == TURB_MODEL::SST) - visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), - solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); + visc_numerics->SetNormal(Normal); + su2double Coord_Reflected[MAXNDIM]; + GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(Point_Normal), + geometry->nodes->GetCoord(iPoint), Coord_Reflected); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), Coord_Reflected); - /*--- Compute and update residual ---*/ + /*--- Primitive variables, and gradient ---*/ - auto residual_v = visc_numerics->ComputeResidual(config); + visc_numerics->SetPrimitive(V_domain, V_inlet); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes->GetGradient_Primitive(iPoint)); - LinSysRes.SubtractBlock(iPoint, residual_v); + /*--- Turbulent kinetic energy ---*/ - /*--- Jacobian contribution for implicit integration ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), + solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual_v.jacobian_i); + /*--- Compute and update residual ---*/ + + auto residual_v = visc_numerics->ComputeResidual(config); + + LinSysRes.SubtractBlock(iPoint, residual_v); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual_v.jacobian_i); + } } END_SU2_OMP_FOR } @@ -2731,6 +2890,9 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, case INC_OUTLET_TYPE::MASS_FLOW_OUTLET: + if (pressure_based) + SU2_MPI::Error("Mass Flow Outlet is currently an unsupported INC_OUTLET_TYPE for pressure based solver.", CURRENT_FUNCTION); + /*--- Retrieve the specified target mass flow at the outlet. ---*/ mDot_Target = config->GetOutlet_Pressure(Marker_Tag)/(config->GetDensity_Ref() * config->GetVelocity_Ref()); @@ -2784,7 +2946,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Beta coefficient from the config file ---*/ - V_outlet[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); + if (!pressure_based) V_outlet[prim_idx.Beta()] = nodes->GetBetaInc2(iPoint); /*--- Cp is needed for Temperature equation. ---*/ @@ -2801,6 +2963,21 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + /*--- Set the edge mass flux ---*/ + + if (pressure_based) { + su2double ProjVelocity = 0.0; + if (dynamic_grid) + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += (V_domain[iDim+prim_idx.Velocity()] - geometry->nodes->GetGridVel(iPoint)[iDim]) * Normal[iDim]; + else + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += V_domain[iDim+prim_idx.Velocity()] * Normal[iDim]; + su2double MeanDensity = 0.5 * (V_domain[prim_idx.Density()] + V_outlet[prim_idx.Density()]); + su2double MassFlux = MeanDensity * ProjVelocity; + conv_numerics->SetMassFlux(MassFlux); + } + /*--- Compute the residual using an upwind scheme ---*/ auto residual = conv_numerics->ComputeResidual(config); @@ -3369,6 +3546,31 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf LoadRestart_impl(geometry, solver, config, val_iter, val_update_geo, Solution, nVar_Restart); + if (pressure_based) { + + /*--- Initialize the edge mass flux array ---*/ + + unsigned long iEdge, iPoint, jPoint; + su2double MeanVelocity[MAXNDIM], MeanDensity; + + for (iEdge = 0; iEdge < geometry[MESH_0]->GetnEdge(); iEdge++) { + + iPoint = geometry[MESH_0]->edges->GetNode(iEdge,0); jPoint = geometry[MESH_0]->edges->GetNode(iEdge,1); + + /*--- Compute average velocities and density between two nodes ---*/ + + for (unsigned short iDim = 0; iDim < nDim; iDim++) + MeanVelocity[iDim] = 0.5 * (nodes->GetVelocity(iPoint, iDim) + nodes->GetVelocity(jPoint, iDim)); + + MeanDensity = 0.5 * (nodes->GetDensity(iPoint) + nodes->GetDensity(jPoint)); + + /*--- Initialize the edge mass flux ---*/ + + EdgeMassFluxes[iEdge] = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + EdgeMassFluxes[iEdge] += MeanDensity * MeanVelocity[iDim] * geometry[MESH_0]->edges->GetNormal(iEdge)[iDim]; + } + } } void CIncEulerSolver::SetFreeStream_Solution(const CConfig *config){ @@ -3408,3 +3610,351 @@ void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, con adj_sol[0] = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop); } } + +void CIncEulerSolver::CorrectPressureGradient(su2double* corrected_grad_pressure, + const su2double* avg_grad_pressure, + const su2double val_pressure_i, + const su2double val_pressure_j, + const su2double* val_edge_vector, + const su2double val_dist_ij_2) { + + /*--- Eq 15.62 F Moukalled, L Mangani M. Darwish OpenFOAM and uFVM book. ---*/ + su2double Proj_Mean_Grad_Pressure_Edge = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Proj_Mean_Grad_Pressure_Edge += avg_grad_pressure[iDim]*val_edge_vector[iDim]; + } + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + corrected_grad_pressure[iDim] = avg_grad_pressure[iDim] - (Proj_Mean_Grad_Pressure_Edge - + (val_pressure_j-val_pressure_i))*val_edge_vector[iDim] / val_dist_ij_2; + } +} + +void CIncEulerSolver::ComputeEdgeMassFluxesRhieChow(CGeometry *geometry, CSolver **solver_container, CConfig *config) { + SU2_ZONE_SCOPED + + /*--- Compute gradients to be used in Rhie Chow interpolation ---*/ + + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetPrimitive_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetPrimitive_Gradient_LS(geometry, config); + } + + unsigned short iDim; + unsigned long iPoint, jPoint; + const su2double *Normal = nullptr, *Coord_i, *Coord_j, *GridVel_i,*GridVel_j; + su2double GradPressure_f[MAXNDIM], GradPressure_avg[MAXNDIM], Edge_Vector[MAXNDIM], dist_ij_2, Coeff_Mom; + + CSolver* poisson_solver = solver_container[POISSON_SOL]; + CVariable* poisson_nodes = poisson_solver->GetNodes(); + + /*--- Mass flux is computed over all edges ---*/ + + for (unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + + iPoint = geometry->edges->GetNode(iEdge,0); jPoint = geometry->edges->GetNode(iEdge,1); + + Normal = geometry->edges->GetNormal(iEdge); + + if (dynamic_grid) { + GridVel_i = geometry->nodes->GetGridVel(iPoint); + GridVel_j = geometry->nodes->GetGridVel(jPoint); + } + + /*--- Correct pressure gradient ---*/ + + Coord_i = geometry->nodes->GetCoord(iPoint); + Coord_j = geometry->nodes->GetCoord(jPoint); + dist_ij_2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; + dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; + } + + /*--- 1. Interpolate the pressure gradient based on node values ---*/ + + for (iDim = 0; iDim < nDim; iDim++) + GradPressure_avg[iDim] = 0.5 * (nodes->GetGradient_Primitive(iPoint,prim_idx.Pressure(),iDim) + nodes->GetGradient_Primitive(jPoint,prim_idx.Pressure(),iDim)); + + /*--- 2. Compute pressure gradient at the face ---*/ + + CorrectPressureGradient(GradPressure_f, GradPressure_avg, nodes->GetPressure(iPoint), nodes->GetPressure(jPoint), Edge_Vector, dist_ij_2); + + /*--- Linearly interpolated coefficient. ---*/ + + Coeff_Mom = 0.5*(poisson_nodes->GetMomCoeff(iPoint) + poisson_nodes->GetMomCoeff(jPoint)); + + /*--- Initialize mass flux ---*/ + + EdgeMassFluxes[iEdge] = 0.0; + + for (iDim = 0; iDim < nDim; iDim++) { + + /*--- Face average mass flux. ---*/ + + su2double meanMassFlux = 0.5 * (nodes->GetDensity(iPoint) * nodes->GetVelocity(iPoint,iDim) + + nodes->GetDensity(jPoint) * nodes->GetVelocity(jPoint,iDim)); + + if (dynamic_grid) { + meanMassFlux -= 0.5 * (nodes->GetDensity(iPoint) * GridVel_i[iDim] + nodes->GetDensity(jPoint) * GridVel_j[iDim]); + } + + /*--- Correction based on Rhie-Chow. ---*/ + + su2double RhieChowCorrection = Coeff_Mom * (GradPressure_f[iDim] - GradPressure_avg[iDim]); + + su2double CorrectedMassFlux = meanMassFlux - RhieChowCorrection; + + /*--- Update edge mass flux ---*/ + + EdgeMassFluxes[iEdge] += CorrectedMassFlux * Normal[iDim]; + + } + } +} + + +void CIncEulerSolver::ApplyPressureVelocityCorrection(CGeometry *geometry, CSolver **solver_container, CConfig *config) { + SU2_ZONE_SCOPED + + /*--- Start of computing the corrections ---*/ + unsigned long iPoint, jPoint, iMarker, iVertex; + unsigned short iDim, KindBC; + su2double Current_Pressure, factor, PCorr_Ref, Vol, delT; + string Marker_Tag; + const su2double *Normal = nullptr; + + bool AutomaticURF = config->GetSIMPLE_Options().AutomaticRelaxationFactors; + + CSolver* poisson_solver = solver_container[POISSON_SOL]; + CVariable* poisson_nodes = poisson_solver->GetNodes(); + + /*--- Combine all pressure corrections into a vector for easy access ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + pressureCorrection[iPoint] = poisson_nodes->GetSolution(iPoint,0); + } + END_SU2_OMP_FOR + + /*--- Define a reference pressure ---*/ + // TODO: look at this, currently copied (but working?) logic from old solver (by Akshay) + unsigned long PRef_Point = 1; + auto Pref_local = geometry->GetGlobal_to_Local_Point(PRef_Point); + PCorr_Ref = 0.0; + if (Pref_local >= 0) + if(geometry->nodes->GetDomain(Pref_local)) + PCorr_Ref = 0.0;//Pressure_Correc[Pref_local]; + + /*--- Compute Velocity Corrections and under relaxation factor for the pressure. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + factor = 0.0; + const auto view = Jacobian.GetBlockView(iPoint, iPoint); + for (iDim = 0; iDim < nDim; iDim++) { + momentumCorrection[iPoint][iDim] = - poisson_nodes->GetMomCoeff(iPoint) * poisson_nodes->GetGradient(iPoint,0,iDim); + if (AutomaticURF) factor += view(iDim, iDim); + } + + if (AutomaticURF) { + Vol = geometry->nodes->GetVolume(iPoint); + delT = nodes->GetDelta_Time(iPoint); + alpha_p[iPoint] = (Vol / delT) / (factor + (Vol / delT)); + } else { + alpha_p[iPoint] = config->GetSIMPLE_Options().Relaxation_Factor_Pressure; + } + + } + END_SU2_OMP_FOR + + // TODO: The HbyA correction is always zero during the first PISO correction, therefore this can be skipped. + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) { + momentumCorrection[iPoint][iDim] += poisson_nodes->GetHbyACorrection(iPoint, iDim); + } + } + END_SU2_OMP_FOR + + /*--- Compute the edge corrections based on the average of the momentum coefficients and the average of the p' gradient. ---*/ + + su2double* Coord_i,* Coord_j; + su2double GradPressure_f[MAXNDIM], GradPressure_avg[MAXNDIM], Edge_Vector[MAXNDIM], dist_ij_2; + for (unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + + iPoint = geometry->edges->GetNode(iEdge,0); jPoint = geometry->edges->GetNode(iEdge,1); + + Normal = geometry->edges->GetNormal(iEdge); + + /*--- Correct pressure deviation (p') gradient ---*/ + + Coord_i = geometry->nodes->GetCoord(iPoint); + Coord_j = geometry->nodes->GetCoord(jPoint); + dist_ij_2 = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; + dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; + } + + /*--- 1. Interpolate the p' gradient based on node values ---*/ + + for (iDim = 0; iDim < nDim; iDim++) + GradPressure_avg[iDim] = 0.5*(poisson_nodes->GetGradient_Primitive(iPoint,0,iDim) + poisson_nodes->GetGradient_Primitive(jPoint,0,iDim)); + + /*--- 2. Compute p' at the face ---*/ + + CorrectPressureGradient(GradPressure_f, GradPressure_avg, poisson_nodes->GetSolution(iPoint, 0), poisson_nodes->GetSolution(jPoint, 0), Edge_Vector, dist_ij_2); + + /*--- Initialize projected velocity and density ---*/ + + su2double ProjMassFluxCorrection = 0.0; + + for (iDim = 0; iDim < nDim; iDim++) { + + su2double MassFluxCorrection = -0.5 * (poisson_nodes->GetMomCoeff(iPoint) + poisson_nodes->GetMomCoeff(jPoint)) * GradPressure_f[iDim]; + + /*--- 2nd piso correction term (HbyA') --- (TODO: this is zero for the first correction and can thus also be skipped) ---*/ + + MassFluxCorrection += 0.5*(poisson_nodes->GetHbyACorrection(iPoint, iDim) + +poisson_nodes->GetHbyACorrection(jPoint, iDim)); + + /*--- Accumulate into the edge mass flux correction ---*/ + + ProjMassFluxCorrection += MassFluxCorrection * Normal[iDim]; + } + + /*--- Set the mass flux correction ---*/ + + EdgeMassFluxCorrection[iEdge] = ProjMassFluxCorrection; + } + + /*--- Reassign strong boundary conditions ---*/ + /*--- For now I only have velocity inlet and fully developed outlet. Will need to add other types of inlet/outlet conditions + * where different treatment of pressure might be needed. Symmetry and Euler wall are weak BCs. ---*/ + for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + KindBC = config->GetMarker_All_KindBC(iMarker); + Marker_Tag = config->GetMarker_All_TagBound(iMarker); + switch (KindBC) { + case EULER_WALL: case SYMMETRY_PLANE: + break; + + /*--- Nothing at MPI boundaries ---*/ + case SEND_RECEIVE: + break; + + /*--- Only a fully developed outlet is implemented. For pressure, a dirichlet + BC has to be applied and no correction is necessary. Velocity has a neumann BC. ---*/ + case OUTLET_FLOW:{ + auto Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag); + switch (Kind_Outlet) { + case INC_OUTLET_TYPE::PRESSURE_OUTLET:{ + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (geometry->nodes->GetDomain(iPoint)) + pressureCorrection[iPoint] = PCorr_Ref; + } + break; + } + //TODO: other outlet types + default: + SU2_MPI::Error("The requested outflow boundary condition has not yet been implemented for the pressure based poisson solver", CURRENT_FUNCTION); + break; + } + break; + } + + /*--- Only a fixed velocity inlet is implemented now. Along with the wall boundaries, + * the velocity is known and thus no correction is necessary.---*/ + case ISOTHERMAL: case HEAT_FLUX: case INLET_FLOW: { + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (geometry->nodes->GetDomain(iPoint)) { + for (iDim = 0; iDim < nDim; iDim++) + momentumCorrection[iPoint][iDim] = 0.0; + alpha_p[iPoint] = 1.0; + } + } + END_SU2_OMP_FOR + break; + } + + /*--- Farfield is treated as a fully developed flow for pressure and a fixed pressure is + * used, thus no correction is necessary. The treatment for velocity depends on whether the + * flow is into the domain or out. If flow is in, a dirichlet bc is applied and no correction + * is made, otherwise a Neumann BC is used and velocity is adjusted. ---*/ + + case FAR_FIELD: + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (geometry->nodes->GetDomain(iPoint)) { + // Check if the boundary condition is an inlet or not + if (nodes->GetStrongBC(iPoint)) { + for (iDim = 0; iDim < nDim; iDim++) + momentumCorrection[iPoint][iDim] = 0.0; + } + pressureCorrection[iPoint] = PCorr_Ref; + } + } + + break; + + default: + SU2_MPI::Error("The requested boundary condition has not yet been implemented for the pressure based poisson solver", CURRENT_FUNCTION); + break; + } + } + + /*--- Apply corrections to the nodal solution ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Pressure corrections ---*/ + + Current_Pressure = nodes->GetPressure(iPoint); + Current_Pressure += alpha_p[iPoint] * (pressureCorrection[iPoint] - PCorr_Ref); + nodes->SetSolution(iPoint, 0, Current_Pressure); + + /*--- Velocity corrections ---*/ + + for (iDim = 0; iDim < nDim; ++iDim) { + nodes->SetSolution(iPoint, iDim + 1, nodes->GetSolution(iPoint,iDim + 1) + momentumCorrection[iPoint][iDim] / nodes->GetDensity(iPoint)); + poisson_nodes->SetMomentumCorrection(iPoint,iDim,momentumCorrection[iPoint][iDim]); + } + + /*--- Update primitive variables ---*/ + + nodes->SetPressure(iPoint); + nodes->SetVelocity(iPoint); + + } + END_SU2_OMP_FOR + + /*--- Add corrections to the edge velocities ---*/ + + for (unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) + EdgeMassFluxes[iEdge] += EdgeMassFluxCorrection[iEdge]; + + /*--- Reset HbyA for next iteration ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPoint; iPoint++) { + for (iDim = 0; iDim < nDim; iDim++) + poisson_nodes->SetHbyACorrection(iPoint, iDim, 0.0); + } + END_SU2_OMP_FOR + + + /*--- periodic communication for both the momentum and the poisson equations as both are now updated ---*/ + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + } + + /*--- Communicate updated velocities and pressure ---*/ + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + +} diff --git a/SU2_CFD/src/solvers/CPoissonSolver.cpp b/SU2_CFD/src/solvers/CPoissonSolver.cpp new file mode 100644 index 00000000000..125b4dd656b --- /dev/null +++ b/SU2_CFD/src/solvers/CPoissonSolver.cpp @@ -0,0 +1,572 @@ +/*! + * \file CPoissonSolver.cpp + * \brief Main subroutines for solving the Poisson equation + * \author F. Palacios, T. Economon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/solvers/CPoissonSolver.hpp" +#include +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/solvers/CScalarSolver.inl" + +/*--- Explicit instantiation of the parent class of CPoissonSolver. ---*/ +template class CScalarSolver; + +CPoissonSolver::CPoissonSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) + : CScalarSolver(geometry, config, false, false, LINEAR_SOLVER_MODE::POISSON) { + SU2_ZONE_SCOPED + + /*--- Dimension of the problem --> pressure deviation is the only conservative variable ---*/ + + nVar = 1; + nPrimVar = 1; + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + /*--- Initialize nVarGrad for deallocation ---*/ + + nVarGrad = nVar; + + /*--- Define geometry constants in the solver structure ---*/ + + nDim = geometry->GetnDim(); + + /*--- Define some structures for locating max residuals ---*/ + + Residual_RMS.resize(nVar,0.0); + Residual_Max.resize(nVar,0.0); + Point_Max.resize(nVar,0); + Point_Max_Coord.resize(nVar,nDim) = su2double(0.0); + + + /*--- Initialization of the structure of the whole Jacobian ---*/ + + if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (poisson equation) MG level: " << iMesh << "." << endl; + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy, false, config->GetKind_Poisson_Linear_Solver_Prec()); + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); + LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + if (ReducerStrategy) EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + + if (config->GetExtraOutput()) { + if (nDim == 2) { nOutputVariables = 13; } + else if (nDim == 3) { nOutputVariables = 19; } + OutputVariables.Initialize(nPoint, nPointDomain, nOutputVariables, 0.0); + OutputHeadingNames = new string[nOutputVariables]; + } + + /*--- Initialize the nodes vector. ---*/ + + nodes = new CPoissonVariable(0.0, nPoint, nDim, nVar, config); + + SetBaseClassPointerToNodes(); + + /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. ---*/ + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + } + + /*--- MPI solution ---*/ + + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + + /*--- Add the solver name. ---*/ + + SolverName = "POISSON"; + +} + + +void CPoissonSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { + SU2_ZONE_SCOPED + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);) + + /*--- Reset pressure corrections to zero for next iteration. ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + nodes->SetSolution(iPoint,0,0.0); + } + END_SU2_OMP_FOR + + /*--- Communicate updated Poisson solution (which should now be zero everywhere) ---*/ + solver_container[POISSON_SOL]->InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + solver_container[POISSON_SOL]->CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + solver_container[POISSON_SOL]->InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + solver_container[POISSON_SOL]->CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + } + + /*--- Compute the gradients only after the solution has been reset to zero ---*/ + CommonPreprocessing(geometry, config, Output); + + /*--- Need to clear EdgeFluxes and Jacobian. ---*/ + if (!Output) { + LinSysRes.SetValZero(); + if (ReducerStrategy) EdgeFluxes.SetValZero(); + Jacobian.SetValZero(); + } + + +} + +void CPoissonSolver::Postprocessing(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh) { + SU2_ZONE_SCOPED + + /*--- Compute gradients of the pressure correction p' so we can use it to find the velocity corrections ---*/ + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) + SetSolution_Gradient_GG(geometry, config,false); + + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) + SetSolution_Gradient_LS(geometry, config,false); + +} + + +void CPoissonSolver::SetMomCoeff(CGeometry *geometry, CSolver **solver_container, CConfig *config, bool periodic, unsigned short iMesh) { + + bool simplec = (config->GetKind_PBIter() == PBITER::SIMPLEC); + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + const CSolver* flow_solution = solver_container[FLOW_SOL]; + const CVariable* flow_nodes = flow_solution->GetNodes(); + + if (implicit) { + + /*--- First sum up the momentum coefficient using the jacobian from given point and it's neighbors. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Self contribution of the coefficient A_p, defined as dR/d(rhou). The jacobian of the momentum + equations is already defined as dR/du so it can be reused. Note that this coefficient should be the same for + all variable directions, therefore just the x-momentum coefficient is taken. ---*/ + + su2double A_p = flow_solution->Jacobian.GetBlockView(iPoint, iPoint)(1,1) / flow_nodes->GetDensity(iPoint); + + /*--- Optionally alter the coefficient using SIMPLEC ---*/ + + su2double Sum_A_nb = 0.0; + + if (simplec) { + for (unsigned long iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { + auto jPoint = geometry->nodes->GetPoint(iPoint,iNeigh); + Sum_A_nb += flow_solution->Jacobian.GetBlockView(iPoint, jPoint)(1,1) / flow_nodes->GetDensity(jPoint); + } + } + + /*--- Add simplec neighbour contributions and optional time dependent term. ---*/ + + su2double Vol = geometry->nodes->GetVolume(iPoint); + su2double delT = flow_nodes->GetDelta_Time(iPoint); + + su2double CorrectedA_p = A_p - Sum_A_nb - config->GetSIMPLE_Options().Transient_Term_Removal_Factor * (Vol / delT); + + /*--- Invert the momentum coefficient to 1/a_p and scale by the volume and density so it can be used as diffusion coefficient in the poisson eq ---*/ + + nodes->SetMomCoeff(iPoint, Vol / CorrectedA_p); + + } + END_SU2_OMP_FOR + } + else { + + SU2_MPI::Error("The definition of the momentum coefficient for an explicit solution is currently only an approximation and is not yet tested.", CURRENT_FUNCTION); + + /* + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + su2double delT = flow_nodes->GetDelta_Time(iPoint); + + su2double Mom_Coeff = delT; + + nodes->SetMomCoeff(iPoint, Mom_Coeff); + } + END_SU2_OMP_FOR + */ + } + + /*--- Insert MPI call here. ---*/ + InitiateComms(geometry, config, MPI_QUANTITIES::MOM_COEFF); + CompleteComms(geometry, config, MPI_QUANTITIES::MOM_COEFF); +} + + +void CPoissonSolver::ComputeHbyA(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { + + unsigned short iDim; + unsigned long iPoint, jPoint, iNeigh; + su2double H, A_p, A_nb; + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + const CSolver* flow_solver = solver_container[FLOW_SOL]; + const CVariable* flow_nodes = flow_solver->GetNodes(); + + /*--- First exchange momentum correction which is required to compute H. ---*/ + InitiateComms(geometry, config, MPI_QUANTITIES::MOM_CORRECTION); + CompleteComms(geometry, config, MPI_QUANTITIES::MOM_CORRECTION); + + if (implicit) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (iDim = 0; iDim < nDim; ++iDim) { + H = 0.0; + A_p = flow_solver->Jacobian.GetBlockView(iPoint, iPoint)(1,1) / flow_nodes->GetDensity(iPoint); + for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { + jPoint = geometry->nodes->GetPoint(iPoint,iNeigh); + A_nb = flow_solver->Jacobian.GetBlockView(iPoint, jPoint)(1,1) / flow_nodes->GetDensity(jPoint); + H -= A_nb * nodes->GetMomentumCorrection(jPoint, iDim); + } + nodes->SetHbyACorrection(iPoint, iDim, H/A_p); + } + } + END_SU2_OMP_FOR + } + else { + SU2_MPI::Error("HbyA is currently not supported for an explicit momentum solver.", CURRENT_FUNCTION); + } + + /*--- Exchange HbyA with MPI call. ---*/ + InitiateComms(geometry, config, MPI_QUANTITIES::HBYA_CORRECTION); + CompleteComms(geometry, config, MPI_QUANTITIES::HBYA_CORRECTION); +} + +void CPoissonSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, + CConfig *config, unsigned short iMesh, unsigned short iRKStep) { + SU2_ZONE_SCOPED + + CNumerics* numerics = numerics_container[VISC_TERM + omp_get_thread_num() * MAX_TERMS]; + + bool pausePreacc = false; + if (ReducerStrategy) + pausePreacc = AD::PausePreaccumulation(); + else + AD::StartNoSharedReading(); + + for (auto color : EdgeColoring) { + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) + for (auto k = 0ul; k < color.size; ++k) { + auto iEdge = color.indices[k]; + Viscous_Residual(iEdge, geometry, solver_container, numerics, config); + } + END_SU2_OMP_FOR + } + + /*--- Restore preaccumulation and adjoint evaluation state. ---*/ + AD::ResumePreaccumulation(pausePreacc); + if (!ReducerStrategy) AD::EndNoSharedReading(); + + if (ReducerStrategy) { + SumEdgeFluxes(geometry); + Jacobian.SetDiagonalAsColumnSum(); + } +} + +void CPoissonSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, + CConfig *config, unsigned short iMesh) { + SU2_ZONE_SCOPED + + su2double *GridVel_i; + + const CSolver* flow_solver = solver_container[FLOW_SOL]; + const CVariable* flow_nodes = flow_solver->GetNodes(); + + const auto& edgeMassFluxes = *(flow_solver->GetEdgeMassFluxes()); + + /*--- flux is computed over all edges ---*/ + + for (auto color : EdgeColoring) { + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) + for (auto k = 0ul; k < color.size; ++k) { + auto iEdge = color.indices[k]; + + auto iPoint = geometry->edges->GetNode(iEdge,0); auto jPoint = geometry->edges->GetNode(iEdge,1); + su2double Normal[MAXNDIM] = {0.0}; + geometry->edges->GetNormal(iEdge, Normal); + + /*--- Add the mass flux to the source term for the poisson equation ---*/ + + auto residual = CNumerics::ResidualType<>(&edgeMassFluxes[iEdge], nullptr, nullptr); + + if (geometry->nodes->GetDomain(iPoint)) LinSysRes.AddBlock(iPoint, residual); + if (geometry->nodes->GetDomain(jPoint)) LinSysRes.SubtractBlock(jPoint, residual); + + /*--- Only for the second pressure correction in the case PISO is used, we need the additional HbyA(u') term ---*/ + // TODO: currently its just set to zero and does not contribute for the first piso correctin but would be nice if this entire block would be skipped otherwise. + su2double MeanHbyA = 0.0; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) + MeanHbyA += 0.5 * (nodes->GetHbyACorrection(iPoint, iDim) + nodes->GetHbyACorrection(jPoint, iDim)) * Normal[iDim]; + + auto residualHbyA = CNumerics::ResidualType<>(&MeanHbyA, nullptr, nullptr); + if (geometry->nodes->GetDomain(iPoint)) LinSysRes.AddBlock(iPoint, residualHbyA); + if (geometry->nodes->GetDomain(jPoint)) LinSysRes.SubtractBlock(jPoint, residualHbyA); + + } + END_SU2_OMP_FOR + } + + /*--- Now add corrections to the previously computed mass fluxes for boundary conditions which alter the mass flux ---*/ + + unsigned short iDim, KindBC; + unsigned long iMarker, iVertex, iPoint; + string Marker_Tag; + su2double MassFlux_corr = 0.0, Normal[MAXNDIM]; + + /*--- Loop boundary edges ---*/ + for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + KindBC = config->GetMarker_All_KindBC(iMarker); + Marker_Tag = config->GetMarker_All_TagBound(iMarker); + + switch (KindBC) { + /*--- Wall boundaries have zero mass flux (irrespective of grid movement) ---*/ + case EULER_WALL: case ISOTHERMAL: case HEAT_FLUX: case SYMMETRY_PLANE: + break; + + /*--- Nothing has to happen at MPI boundaries*/ + case SEND_RECEIVE: + break; + + case INLET_FLOW: + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (!geometry->nodes->GetDomain(iPoint)) continue; + + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + + MassFlux_corr = 0.0; + if (dynamic_grid) { + GridVel_i = geometry->nodes->GetGridVel(iPoint); + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * (flow_nodes->GetVelocity(iPoint, iDim) - GridVel_i[iDim]) * Normal[iDim]; + } + else + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * flow_nodes->GetVelocity(iPoint, iDim) * Normal[iDim]; + + auto residual = CNumerics::ResidualType<>(&MassFlux_corr, nullptr, nullptr); + + if (geometry->nodes->GetDomain(iPoint)) LinSysRes.AddBlock(iPoint, residual); + + } + break; + + case FAR_FIELD: + /*--- Treat the farfield as a fully developed outlet for pressure. ---*/ + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (geometry->nodes->GetDomain(iPoint)) { + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + + if (dynamic_grid) + GridVel_i = geometry->nodes->GetGridVel(iPoint); + + MassFlux_corr = 0.0; + if (dynamic_grid) + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * (flow_nodes->GetVelocity(iPoint, iDim) - GridVel_i[iDim]) * Normal[iDim]; + else + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * flow_nodes->GetVelocity(iPoint, iDim) * Normal[iDim]; + + auto residual = CNumerics::ResidualType<>(&MassFlux_corr, nullptr, nullptr); + LinSysRes.AddBlock(iPoint, residual); + + } + } + break; + + + case OUTLET_FLOW:{ + /*--- Note I am assuming a fully developed outlet, thus the pressure value is prescribed + * -- and a dirichlet bc has to be applied along outlet faces. The Massflux, which forms the RHS + * -- of the equation, is set to zero to enforce the dirichlet bc. ---*/ + + auto Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag); + + switch (Kind_Outlet) { + case INC_OUTLET_TYPE::PRESSURE_OUTLET: + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (geometry->nodes->GetDomain(iPoint)) { + geometry->vertex[iMarker][iVertex]->GetNormal(Normal); + + if (dynamic_grid) + GridVel_i = geometry->nodes->GetGridVel(iPoint); + + MassFlux_corr = 0.0; + if (dynamic_grid) + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * (flow_nodes->GetVelocity(iPoint, iDim) - GridVel_i[iDim]) * Normal[iDim]; + else + for (iDim = 0; iDim < nDim; iDim++) + MassFlux_corr -= flow_nodes->GetDensity(iPoint) * (flow_nodes->GetVelocity(iPoint, iDim)) * Normal[iDim]; + + auto residual = CNumerics::ResidualType<>(&MassFlux_corr, nullptr, nullptr); + + if (geometry->nodes->GetDomain(iPoint)) LinSysRes.AddBlock(iPoint, residual); + } + } + break; + default: + SU2_MPI::Error("Requested type of outlet boundary condition not available", CURRENT_FUNCTION); + break; + } + break; + } + + default: + SU2_MPI::Error("Invalid boundary condition for flux correction", CURRENT_FUNCTION); + break; + + } + } + +} + +void CPoissonSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { + SU2_ZONE_SCOPED + + /*--- No actual time integration is done here. The routine is used as a means to solve the linear equation + * resulting from the poisson equation. The linear system is solved using the jacobian matrix in a way + * consistent with the rest of the code. The time step is set to zero and no under-relaxation is applied to the + * jacobian matrix. ---*/ + + /*--- Local residual variables for current thread ---*/ + su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; + unsigned long idxMax[MAXNVAR] = {0}; + + SetResToZero(); + + /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/ + SU2_OMP_FOR_(schedule(static,omp_chunk_size) SU2_NOWAIT) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Multigrid contribution to residual. ---*/ + su2double *local_Res_TruncError = nodes->GetResTruncError(iPoint); + + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + LinSysRes(iPoint, iVar) = - (LinSysRes(iPoint, iVar) + local_Res_TruncError[iVar] ); + LinSysSol(iPoint, iVar) = 0.0; + + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint, iVar, LinSysRes(iPoint, iVar), resRMS, resMax, idxMax); + } + } + END_SU2_OMP_FOR + + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); + + /*--- Solve or smooth the linear system. ---*/ + + SU2_OMP_FOR_(schedule(static,OMP_MIN_SIZE) SU2_NOWAIT) + for (unsigned long iPoint = nPointDomain; iPoint < nPoint; iPoint++) { + LinSysRes.SetBlock_Zero(iPoint); + LinSysSol.SetBlock_Zero(iPoint); + } + END_SU2_OMP_FOR + + auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + SetIterLinSolver(iter); + SetResLinSolver(System.GetResidual()); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + nodes->AddSolution(iPoint, iVar, LinSysSol(iPoint,iVar)); + } + } + END_SU2_OMP_FOR + + InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + +} + +void CPoissonSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + /*--- Zero flux (Neumann) BC on pressure ---*/ +} + +void CPoissonSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + + unsigned long iVertex, iPoint; + su2double pressureDeviation = 0.0; + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; + /*--- The farfield boundary is considered as an inlet-outlet boundary, where flow + * can either enter or leave. For pressure, it is treated as a fully developed flow + * and a dirichlet BC is applied. For velocity, based on the sign of massflux, either + * a dirichlet or a neumann BC is applied (in correction routine). ---*/ + + LinSysRes.SetBlock_Zero(iPoint); + + nodes->SetSolution(iPoint, &pressureDeviation); + nodes->SetSolution_Old(iPoint,&pressureDeviation); + Jacobian.DeleteValsRowi(iPoint, 0); + } +} + +void CPoissonSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + /*--- Zero flux (Neumann) BC on pressure ---*/ +} + +void CPoissonSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + unsigned long iVertex, iPoint; + su2double pressureDeviation = 0.0; + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; + + /*--- apply a dirichlet boundary condition as pressure is prescribed*/ + + LinSysRes.SetBlock_Zero(iPoint); + + nodes->SetSolution(iPoint, &pressureDeviation); + nodes->SetSolution_Old(iPoint,&pressureDeviation); + Jacobian.DeleteValsRowi(iPoint, 0); + } +} \ No newline at end of file diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 5662aa46ed4..ac90a1ad585 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1399,6 +1399,18 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE::DOUBLE; break; + case MPI_QUANTITIES::MOM_COEFF: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE::DOUBLE; + break; + case MPI_QUANTITIES::MOM_CORRECTION: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE::DOUBLE; + break; + case MPI_QUANTITIES::HBYA_CORRECTION: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE::DOUBLE; + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); @@ -1556,6 +1568,17 @@ void CSolver::InitiateComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_time_n1(iPoint, iVar); break; + case MPI_QUANTITIES::MOM_COEFF: + bufDSend[buf_offset] = base_nodes->GetMomCoeff(iPoint); + break; + case MPI_QUANTITIES::MOM_CORRECTION: + for (iDim = 0; iDim < nDim; iDim++) + bufDSend[buf_offset+iDim] = base_nodes->GetMomentumCorrection(iPoint, iDim); + break; + case MPI_QUANTITIES::HBYA_CORRECTION: + for (iDim = 0; iDim < nDim; iDim++) + bufDSend[buf_offset+iDim] = base_nodes->GetHbyACorrection(iPoint, iDim); + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); @@ -1712,6 +1735,17 @@ void CSolver::CompleteComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) base_nodes->Set_Solution_time_n1(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; + case MPI_QUANTITIES::MOM_COEFF: + base_nodes->SetMomCoeff(iPoint, bufDRecv[buf_offset]); + break; + case MPI_QUANTITIES::MOM_CORRECTION: + for (iDim = 0; iDim < nDim; iDim++) + base_nodes->SetMomentumCorrection(iPoint, iDim, bufDRecv[buf_offset+iDim]); + break; + case MPI_QUANTITIES::HBYA_CORRECTION: + for (iDim = 0; iDim < nDim; iDim++) + base_nodes->SetHbyACorrection(iPoint, iDim, bufDRecv[buf_offset+iDim]); + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index ee798c1384b..f850e13bd18 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -31,6 +31,7 @@ #include "../../include/solvers/CIncEulerSolver.hpp" #include "../../include/solvers/CNSSolver.hpp" #include "../../include/solvers/CIncNSSolver.hpp" +#include "../../include/solvers/CPoissonSolver.hpp" #include "../../include/solvers/CNEMOEulerSolver.hpp" #include "../../include/solvers/CNEMONSSolver.hpp" #include "../../include/solvers/CTurbSASolver.hpp" @@ -69,6 +70,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon case MAIN_SOLVER::INC_EULER: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_EULER, solver, geometry, config, iMGLevel); solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::EULER: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::EULER, solver, geometry, config, iMGLevel); @@ -81,6 +83,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::NAVIER_STOKES: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel); @@ -103,6 +106,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon solver[TURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel); solver[TRANS_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::TRANSITION, solver, geometry, config, iMGLevel); solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::HEAT_EQUATION: solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); @@ -144,6 +148,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon solver[ADJFLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FLOW, solver, geometry, config, iMGLevel); solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); solver[ADJRAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); @@ -154,6 +159,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon solver[ADJRAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel); solver[SPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::SPECIES, solver, geometry, config, iMGLevel); solver[ADJSPECIES_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_SPECIES, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::DISC_ADJ_INC_RANS: solver[FLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel); @@ -166,6 +172,7 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon solver[ADJTURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_TURB, solver, geometry, config, iMGLevel); solver[RAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel); solver[ADJRAD_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel); + solver[POISSON_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::POISSON, solver, geometry, config, iMGLevel); break; case MAIN_SOLVER::DISC_ADJ_HEAT: solver[HEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel); @@ -331,6 +338,11 @@ CSolver* CSolverFactory::CreateSubSolver(SUB_SOLVER_TYPE kindSolver, CSolver **s } metaData.integrationType = INTEGRATION_TYPE::DEFAULT; break; + case SUB_SOLVER_TYPE::POISSON: + if (config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED) + genericSolver = new CPoissonSolver(geometry, config, iMGLevel); + metaData.integrationType = INTEGRATION_TYPE::SINGLEGRID; + break; default: SU2_MPI::Error("No proper allocation found for requested sub solver", CURRENT_FUNCTION); break; diff --git a/SU2_CFD/src/variables/CFlowVariable.cpp b/SU2_CFD/src/variables/CFlowVariable.cpp index 7828b1afbb0..c9f3524ba05 100644 --- a/SU2_CFD/src/variables/CFlowVariable.cpp +++ b/SU2_CFD/src/variables/CFlowVariable.cpp @@ -50,7 +50,7 @@ CFlowVariable::CFlowVariable(unsigned long npoint, unsigned long ndim, unsigned Primitive.resize(nPoint, nPrimVar) = su2double(0.0); - if (config->GetMUSCL_Flow() || config->GetViscous() || config->GetContinuous_Adjoint()) { + if (config->GetMUSCL_Flow() || config->GetViscous() || config->GetContinuous_Adjoint() || config->GetKind_Incomp_System() == INCOMP_SYSTEM::PRESSURE_BASED) { Gradient_Primitive.resize(nPoint, nPrimVarGrad, nDim, 0.0); } diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 5274df76644..48bb743ebbf 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -71,6 +71,10 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci if (config->GetStreamwise_Periodic_Temperature()) Streamwise_Periodic_RecoveredTemperature.resize(nPoint) = su2double(0.0); } + + /*--- Allocate strong BC vector for pressure-based solver ---*/ + + strongBC.resize(nPoint); } bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { diff --git a/SU2_CFD/src/variables/CPoissonVariable.cpp b/SU2_CFD/src/variables/CPoissonVariable.cpp new file mode 100644 index 00000000000..981768aeaed --- /dev/null +++ b/SU2_CFD/src/variables/CPoissonVariable.cpp @@ -0,0 +1,57 @@ +/*! + * \file CPoissonVariable.cpp + * \brief Definition of the variables for poisson equation problems. + * \author F. Palacios, T. Economon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/variables/CPoissonVariable.hpp" + +CPoissonVariable::CPoissonVariable(su2double value, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) + : CScalarVariable(npoint, ndim, nvar, config) { + + /*--- Initialization ---*/ + + Solution = value; + Solution_Old = value; + + /*--- Allocate residual structures ---*/ + + Res_TruncError.resize(nPoint, nVar) = su2double(0.0); + + /*--- Only for residual smoothing (multigrid) ---*/ + + for (unsigned long iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) { + if (config->GetMGOptions().MG_CorrecSmooth[iMesh] > 0) { + Residual_Sum.resize(nPoint, nVar); + Residual_Old.resize(nPoint, nVar); + break; + } + } + + /*--- Initialize momentum coefficient and HbyA ---*/ + MomCoeff.resize(nPoint) = su2double(0.0); + MomentumCorrection.resize(nPoint, nDim) = su2double(0.0); + HbyACorrection.resize(nPoint, nDim) = su2double(0.0); + +} \ No newline at end of file diff --git a/TestCases/incomp_euler/naca0012/incomp_pb_NACA0012.cfg b/TestCases/incomp_euler/naca0012/incomp_pb_NACA0012.cfg new file mode 100644 index 00000000000..48dc770a493 --- /dev/null +++ b/TestCases/incomp_euler/naca0012/incomp_pb_NACA0012.cfg @@ -0,0 +1,109 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Incompressible flow hydrofoil 5 degrees % +% Author: Francisco Palacios % +% Institution: Stanford University % +% Date: 09/18/2011 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_EULER +KIND_INCOMP_SYSTEM= PRESSURE_BASED +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_INIT= 998.2 +INC_VELOCITY_INIT= ( 1.775, 0.0, 0.0 ) +INC_INLET_TYPE= VELOCITY_INLET +INC_INLET_DAMPING= 0.1 +INC_OUTLET_TYPE= PRESSURE_OUTLET +INC_OUTLET_DAMPING= 0.1 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 + +% ----------------------- BOUNDARY CONDITION DEFINITION -----------------------% +% +MARKER_EULER= ( airfoil, lower_wall, upper_wall ) +MARKER_INLET= ( inlet, 0.0, 1.775, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 0.0 ) +MARKER_PLOTTING= ( airfoil ) +MARKER_MONITORING= ( airfoil ) + +% ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 9999 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 25 + +% ----------------------- PRESSURE BASED PARAMETERS ---------------------------% +% +RELAXATION_FACTOR_PRESSURE= 1.0 +TRANSIENT_TERM_REMOVAL_FACTOR= 0.0 +KIND_PB_ITER= SIMPLE +PISO_CORRECTIONS = 2 +POISSON_LINEAR_SOLVER= FGMRES +POISSON_LINEAR_SOLVER_PREC= LU_SGS +POISSON_LINEAR_SOLVER_ERROR= 1E-6 +POISSON_LINEAR_SOLVER_ITER= 1000 + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 0 +MGCYCLE= W_CYCLE +MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) +MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= UDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +VENKAT_LIMITER_COEFF= 0.0002 +JST_SENSOR_COEFF= ( 0.5, 0.04 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -10 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 50 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_NACA0012_5deg_6814.su2 +MESH_FORMAT= SU2 +SOLUTION_FILENAME= solution_flow +TABULAR_FORMAT= CSV +CONV_FILENAME= history_pb +RESTART_FILENAME= restart_flow +VOLUME_FILENAME= flow_pb +SURFACE_FILENAME= surface_flow_pb +OUTPUT_WRT_FREQ= 100 +SCREEN_OUTPUT= (INNER_ITER, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_MASSFLUX, LIFT, DRAG) diff --git a/TestCases/incomp_navierstokes/cylinder/incomp_pb_cylinder.cfg b/TestCases/incomp_navierstokes/cylinder/incomp_pb_cylinder.cfg new file mode 100644 index 00000000000..259c83828f4 --- /dev/null +++ b/TestCases/incomp_navierstokes/cylinder/incomp_pb_cylinder.cfg @@ -0,0 +1,102 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Steady incompressible laminar flow around a cylinder % +% Author: Francisco Palacios % +% Institution: Stanford University % +% Date: 2012.03.14 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_INCOMP_SYSTEM= PRESSURE_BASED +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_INIT= 998.2 +INC_VELOCITY_INIT= ( 0.000008, 0.0, 0.0 ) + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.798E-3 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( cylinder, 0.0 ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING= ( cylinder ) +MARKER_MONITORING= ( cylinder ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 50.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 5000 +VENKAT_LIMITER_COEFF= 0.01 + +% ----------------------- PRESSURE BASED PARAMETERS ---------------------------% +% +RELAXATION_FACTOR_PRESSURE= 1.0 +TRANSIENT_TERM_REMOVAL_FACTOR= 0.0 +KIND_PB_ITER= SIMPLE +PISO_CORRECTIONS = 4 +POISSON_LINEAR_SOLVER= FGMRES +POISSON_LINEAR_SOLVER_PREC= LU_SGS +POISSON_LINEAR_SOLVER_ERROR= 1E-6 +POISSON_LINEAR_SOLVER_ITER= 1000 + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 0 +MGCYCLE= V_CYCLE +MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) +MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= UDS +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW= NONE +JST_SENSOR_COEFF= ( 0.5, 0.04 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -10 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_cylinder_lam.su2 +MESH_FORMAT= SU2 +SOLUTION_FILENAME= solution_flow +TABULAR_FORMAT= CSV +CONV_FILENAME= history_pb +RESTART_FILENAME= restart_flow +VOLUME_FILENAME= flow_pb +SURFACE_FILENAME= surface_flow +OUTPUT_WRT_FREQ= 10 +SCREEN_OUTPUT= (INNER_ITER, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_MASSFLUX, LIFT, DRAG) diff --git a/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_liddrivencavity.cfg b/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_liddrivencavity.cfg new file mode 100644 index 00000000000..eb5b9c17bb6 --- /dev/null +++ b/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_liddrivencavity.cfg @@ -0,0 +1,106 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Lid driven flow inside a cavity % +% Author: T. Aalbers % +% Date: 2025.06.25 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1.0 +INC_VELOCITY_INIT= ( 0.0, 0.0, 0.0 ) + +INC_DENSITY_REF= 1.0 +INC_VELOCITY_REF= 1.0 +INC_NONDIM= REFERENCE_VALUES + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 2.5E-3 + +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% + +% Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, +% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, EXTERNAL) +SURFACE_MOVEMENT= MOVING_WALL +MARKER_MOVING= Top +SURFACE_MOTION_ORIGIN = 0.0 0.0 0.0 +SURFACE_ROTATION_RATE = 0.0 0.0 0.0 +SURFACE_TRANSLATION_RATE = 1.0 0.0 0.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_AREA= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( Top, 0.0, Bottom, 0.0 ) +MARKER_ISOTHERMAL= ( Left, 461.04, Right, 115.26 ) +MARKER_PLOTTING= ( Top, Left, Right, Bottom ) +MARKER_MONITORING= ( Left, Right, Top ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 1000 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 0.95, 1.01, 10, 1000, 0.001, 0) +MAX_DELTA_TIME= 1E6 +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 99999 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 1E-06 +LINEAR_SOLVER_ITER= 25 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -10 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= square_50x50.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= solution_flow +SOLUTION_ADJ_FILENAME= solution_adj +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow +RESTART_ADJ_FILENAME= restart_adj +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 100 +SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y) diff --git a/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_pb_liddrivencavity.cfg b/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_pb_liddrivencavity.cfg new file mode 100644 index 00000000000..0276b922888 --- /dev/null +++ b/TestCases/incomp_navierstokes/lid_driven_cavity/incomp_pb_liddrivencavity.cfg @@ -0,0 +1,114 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Lid driven flow inside a cavity % +% Author: T. Aalbers % +% Date: 2025.06.25 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_INCOMP_SYSTEM = PRESSURE_BASED +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1.0 +INC_VELOCITY_INIT= ( 0.0, 0.0, 0.0 ) + +INC_DENSITY_REF= 1.0 +INC_VELOCITY_REF= 1.0 +INC_NONDIM= REFERENCE_VALUES + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 2.5E-3 + +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% + +% Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, +% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, EXTERNAL) +SURFACE_MOVEMENT= MOVING_WALL +MARKER_MOVING= Top +SURFACE_MOTION_ORIGIN = 0.0 0.0 0.0 +SURFACE_ROTATION_RATE = 0.0 0.0 0.0 +SURFACE_TRANSLATION_RATE = 1.0 0.0 0.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_AREA= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( Top, 0.0, Bottom, 0.0, Left, 0.0, Right, 0.0) +MARKER_PLOTTING= ( Top, Left, Right, Bottom ) +MARKER_MONITORING= ( Left, Right, Top ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 100 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 0.95, 1.01, 10, 1000, 0.001, 0) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 99999 + +% ----------------------- PRESSURE BASED PARAMETERS ---------------------------% +% +RELAXATION_FACTOR_PRESSURE= 1.0 +TRANSIENT_TERM_REMOVAL_FACTOR= 0.0 +KIND_PB_ITER= SIMPLE +PISO_CORRECTIONS = 2 +POISSON_LINEAR_SOLVER= FGMRES +POISSON_LINEAR_SOLVER_PREC= ILU +POISSON_LINEAR_SOLVER_ERROR= 1E-6 +POISSON_LINEAR_SOLVER_ITER= 1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= LU_SGS +LINEAR_SOLVER_ERROR= 1E-06 +LINEAR_SOLVER_ITER= 1000 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= UDS +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -10 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= square_50x50.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= solution_flow +SOLUTION_ADJ_FILENAME= solution_adj +TABULAR_FORMAT= CSV +CONV_FILENAME= history_pb +RESTART_FILENAME= restart_flow +RESTART_ADJ_FILENAME= restart_adj +VOLUME_FILENAME= flow_pb +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 20 +SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, LINSOL) \ No newline at end of file diff --git a/TestCases/incomp_rans/rough_flatplate/pb_rough_flatplate_incomp.cfg b/TestCases/incomp_rans/rough_flatplate/pb_rough_flatplate_incomp.cfg new file mode 100644 index 00000000000..fd44a2a501d --- /dev/null +++ b/TestCases/incomp_rans/rough_flatplate/pb_rough_flatplate_incomp.cfg @@ -0,0 +1,130 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Turbulent flow over rough flat plate with zero % +% pressure gradient % +% Author: Akshay Koodly % +% Date: 2020.07.07 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +KIND_INCOMP_SYSTEM= PRESSURE_BASED +KIND_TURB_MODEL= SA +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_ENERGY_EQUATION = NO +INC_DENSITY_INIT= 1.32905 +INC_VELOCITY_INIT= ( 69.4448, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +INC_NONDIM= INITIAL_VALUES +INC_DENSITY_REF= 1.0 +INC_VELOCITY_REF= 1.0 +INC_TEMPERATURE_REF = 1.0 +INC_INLET_TYPE= VELOCITY_INLET + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.84592e-05 +MU_REF= 1.716E-5 +MU_T_REF= 273.15 +SUTHERLAND_CONSTANT= 110.4 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 2.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall, 0.0 ) +WALL_ROUGHNESS = (wall, 0.000246) +MARKER_INLET= ( inlet, 300.0, 69.4448, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 0.0, farfield, 0.0 ) +INC_OUTLET_TYPE= PRESSURE_OUTLET,PRESSURE_OUTLET +MARKER_SYM= ( symmetry ) +MARKER_PLOTTING= ( wall ) +MARKER_MONITORING= ( wall ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.1, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 1000 + +% ----------------------- PRESSURE BASED PARAMETERS ---------------------------% +% +RELAXATION_FACTOR_PRESSURE= 1.0 +TRANSIENT_TERM_REMOVAL_FACTOR= 0.0 +KIND_PB_ITER= SIMPLE +PISO_CORRECTIONS = 2 +POISSON_LINEAR_SOLVER= FGMRES +POISSON_LINEAR_SOLVER_PREC= ILU +POISSON_LINEAR_SOLVER_ERROR= 1E-6 +POISSON_LINEAR_SOLVER_ITER= 1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 1E-12 +LINEAR_SOLVER_ITER= 20 + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.1 +ADJ_SHARP_LIMITER_COEFF= 3.0 +REF_SHARP_EDGES= 3.0 +SENS_REMOVE_SHARP= NO + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= UDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +TIME_DISCRE_TURB= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_VELOCITY-X +CONV_RESIDUAL_MINVAL= -14 +CONV_STARTITER= 10 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_flatplate_turb_137x97.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= restart_flow_pb +TABULAR_FORMAT= CSV +OUTPUT_FILES= RESTART, PARAVIEW, SURFACE_PARAVIEW +CONV_FILENAME= history +RESTART_FILENAME= restart_flow_pb +VOLUME_FILENAME= flow_pb +SURFACE_FILENAME= surface_flow_pb +OUTPUT_WRT_FREQ= 100, 50, 50 +SCREEN_OUTPUT= (WALL_TIME,INNER_ITER, RMS_VELOCITY-X, RMS_NU_TILDE, LIFT,DRAG) +WRT_FORCES_BREAKDOWN= YES diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ff9b8c66706..43535cb52d7 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -650,6 +650,14 @@ def main(): inc_euler_naca0012.test_vals = [-6.067964, -5.125607, 0.525745, 0.008772] test_list.append(inc_euler_naca0012) + # NACA0012 Hydrofoil + inc_euler_naca0012_pb = TestCase('inc_euler_naca0012_pb') + inc_euler_naca0012_pb.cfg_dir = "incomp_euler/naca0012" + inc_euler_naca0012_pb.cfg_file = "incomp_pb_NACA0012.cfg" + inc_euler_naca0012_pb.test_iter = 20 + inc_euler_naca0012_pb.test_vals = [-4.154445, 4.919399, 0.468631, 0.010040] + test_list.append(inc_euler_naca0012_pb) + # C-D nozzle with pressure inlet and mass flow outlet inc_nozzle = TestCase('inc_nozzle') inc_nozzle.cfg_dir = "incomp_euler/nozzle" @@ -678,6 +686,14 @@ def main(): inc_lam_cylinder.test_vals = [-4.156113, -3.553508, -0.024563, 5.105605] test_list.append(inc_lam_cylinder) + # Laminar cylinder, pressure-based + inc_lam_cylinder_pb = TestCase('inc_lam_cylinder_pb') + inc_lam_cylinder_pb.cfg_dir = "incomp_navierstokes/cylinder" + inc_lam_cylinder_pb.cfg_file = "incomp_pb_cylinder.cfg" + inc_lam_cylinder_pb.test_iter = 10 + inc_lam_cylinder_pb.test_vals = [-3.484504, -3.775524, 0.012091, 6.190144] + test_list.append(inc_lam_cylinder_pb) + # Laminar sphere, Re=1. Last column: Cd=24/Re inc_lam_sphere = TestCase('inc_lam_sphere') inc_lam_sphere.cfg_dir = "incomp_navierstokes/sphere" diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index aa22e30d8ae..d2eedb4f054 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -402,6 +402,14 @@ def main(): inc_euler_naca0012.test_vals = [-5.988713, -5.020635, 0.522968, 0.008854] test_list.append(inc_euler_naca0012) + # NACA0012 Hydrofoil, pressure-based + inc_euler_naca0012_pb = TestCase('inc_euler_naca0012_pb') + inc_euler_naca0012_pb.cfg_dir = "incomp_euler/naca0012" + inc_euler_naca0012_pb.cfg_file = "incomp_pb_NACA0012.cfg" + inc_euler_naca0012_pb.test_iter = 20 + inc_euler_naca0012_pb.test_vals = [-4.154445, -4.919400, 0.468630, 0.010040] + test_list.append(inc_euler_naca0012_pb) + # C-D nozzle with pressure inlet and mass flow outlet inc_nozzle = TestCase('inc_nozzle') inc_nozzle.cfg_dir = "incomp_euler/nozzle" @@ -429,6 +437,14 @@ def main(): inc_lam_cylinder.test_vals = [-4.161215, -3.573002, 0.019888, 4.945923] test_list.append(inc_lam_cylinder) + # Laminar cylinder, pressure-based + inc_lam_cylinder_pb = TestCase('inc_lam_cylinder_pb') + inc_lam_cylinder_pb.cfg_dir = "incomp_navierstokes/cylinder" + inc_lam_cylinder_pb.cfg_file = "incomp_pb_cylinder.cfg" + inc_lam_cylinder_pb.test_iter = 10 + inc_lam_cylinder_pb.test_vals = [-3.484442, -3.775560, 0.012341, 6.190300] + test_list.append(inc_lam_cylinder_pb) + # Buoyancy-driven cavity inc_buoyancy = TestCase('inc_buoyancy') inc_buoyancy.cfg_dir = "incomp_navierstokes/buoyancy_cavity" diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index aefbc66f815..8517f5bd6fb 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -128,6 +128,14 @@ def main(): von_karman_cylinder.test_vals = [-7.845765, -7.681042, -8.736704, -0.002581, 1.423652] test_list.append(von_karman_cylinder) + # Lid Driven Cavity Flow (Re=400), pressure-based + lid_driven_cavity = TestCase('lid_driven_cavity') + lid_driven_cavity.cfg_dir = "../Tutorials/incompressible_flow/Inc_Lid_Driven_Cavity" + lid_driven_cavity.cfg_file = "incomp_pb_liddrivencavity.cfg" + lid_driven_cavity.test_iter = 2 + lid_driven_cavity.test_vals = [-5.798956, -4.393539, -5.070494] + test_list.append(lid_driven_cavity) + ### Species Transport diff --git a/config_template.cfg b/config_template.cfg index 535ae564510..82600b640c9 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -347,7 +347,42 @@ INC_OUTLET_DAMPING= 0.1 BULK_MODULUS= 1.42E5 % Epsilon^2 multipier in Beta calculation for incompressible preconditioner. BETA_FACTOR= 4.1 + +% -------------- PRESSURE-BASED INCOMPRESSIBLE SOLVER DEFINITION ---------------% +% +% Incompressible solver. Options are: DENSITY_BASED (default), PRESSURE_BASED +KIND_INCOMP_SYSTEM= PRESSURE_BASED +% +% Pressure based method. Options are: SIMPLE (default), SIMPLEC +KIND_PB_ITER= SIMPLE +% +% The relaxation factor for the pressure correction (default = 1.0) +RELAXATION_FACTOR_PRESSURE= 1.0 +% +% Decide if automatic relaxation factors for pressure and momentum should be used +% or not. If set to NO, the value from RELAXATION_FACTOR_PRESSURE is used. default NO +USE_AUTOMATIC_RELAXATION_FACTORS= YES +% +% Heuristic parameter influencing the effect of the pseudo time derivative in the +% momentum coefficients for the pressure correction equation (float [0,1], default=0.0) +TRANSIENT_TERM_REMOVAL_FACTOR= 0.0 % +% Number of pressure corrections in the SIMPLE algorithm. >1 leads to the PISO +% algorithm. (default = 1 = SIMPLE) +PISO_CORRECTIONS = 2 +% +% Linear solver for poisson equation (same options as LINEAR_SOLVER) +POISSON_LINEAR_SOLVER= FGMRES +% +% Preconditioner of the poisson linear solver (same options as LINEAR_SOLVER_PREC) +POISSON_LINEAR_SOLVER_PREC= ILU +% +% Minimum error of the linear solver for poisson eq. +POISSON_LINEAR_SOLVER_ERROR= 1E-6 +% +% Max. Number of iterations for the poisson linear solver +POISSON_LINEAR_SOLVER_ITER= 1000 + % ----------------------------- SOLID ZONE HEAT VARIABLES-----------------------% % % Thermal conductivity used for heat equation