From 64d42135a46883cf33ccbb866fb82eb6a33269ab Mon Sep 17 00:00:00 2001 From: Thijs Aalbers Date: Thu, 7 May 2026 14:42:43 +0200 Subject: [PATCH 1/4] Initial configuration options and a test case --- Common/include/CConfig.hpp | 63 +++- Common/include/option_structure.hpp | 42 ++- Common/src/CConfig.cpp | 21 ++ .../lid_driven_cavity/lid_driven_cavity.cfg | 272 ++++++++++++++++++ 4 files changed, 395 insertions(+), 3 deletions(-) create mode 100644 TestCases/incomp_pressure_based/lid_driven_cavity/lid_driven_cavity.cfg diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 71f780bbed0..9eff410c418 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. */ @@ -538,6 +539,7 @@ class CConfig { Kind_TimeIntScheme_AdjTurb, /*!< \brief Time integration for the adjoint turbulence model. */ Kind_TimeIntScheme_Species, /*!< \brief Time integration for the species model. */ Kind_TimeIntScheme_Heat, /*!< \brief Time integration for the wave equations. */ + Kind_TimeIntScheme_Poisson, /*!< \brief Time integration for the poisson (pressure correction) equation/ */ Kind_TimeStep_Heat, /*!< \brief Time stepping method for the (fvm) heat equation. */ n_Datadriven_files; @@ -587,6 +589,11 @@ class CConfig { Kind_Upwind_Heat, /*!< \brief Upwind scheme for the heat transfer model. */ Kind_Upwind_Template; /*!< \brief Upwind scheme for the template model. */ + ENUM_PBITER /*< \brief Pressure-based solver for incompressible flows. */ + Kind_PBIter; + ENUM_INCOMP_SYSTEM + Kind_Incomp_System; + 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.*/ @@ -640,6 +647,7 @@ class CConfig { su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ 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. */ @@ -649,6 +657,8 @@ 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. */ + su2double RCFactor; /*!< \brief Relaxation for the Rhie-Chow interpolation contribution. */ + su2double Relaxation_Factor_PBFlow; /*!< \brief Relaxation coefficient of the flow corrections in the PB solver. */ 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. */ @@ -3992,6 +4002,30 @@ class CConfig { */ ENUM_REGIME GetKind_Regime(void) const { return Kind_Regime; } + /*! + * \brief Kind of incompressible solver formulation. + * \return Kind of incompressible solver. + */ + ENUM_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. + */ + ENUM_PBITER GetKind_PBIter(void) const { return Kind_PBIter; } + + /*! + * \brief Set the kind of incompressible solver formulation that is used. + * \param[in] val_system - the type of system to use. + */ + void SetIncomp_System(unsigned short val_system); + + /*! + * \brief Set the pressure based iteration method. + * \param[in] val_PBIter - The iteration method for the pressure based solver. + */ + void SetPBIter(unsigned short val_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. @@ -4366,6 +4400,12 @@ class CConfig { */ 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. @@ -4399,6 +4439,18 @@ class CConfig { * \return Relaxation factor. */ su2double GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; } + + /*! + * \brief Get the relaxation coefficient of the flow correction for PB solver. + * \return relaxation coefficient of the flow correction for PB solver + */ + su2double GetRelaxation_Factor_PBFlow(void) const { return Relaxation_Factor_PBFlow; } + + /*! + * \brief Get the relaxation coefficient for the Rhie-Chow interpolation in the PB solver. + * \return relaxation coefficient of the Rhie-Chow interpolation. + */ + su2double GetRCFactor(void) const { return RCFactor; } /*! * \brief Get the relaxation factor for solution updates of adjoint solvers. @@ -4655,6 +4707,15 @@ class CConfig { */ unsigned short GetKind_TimeIntScheme(void) const { return Kind_TimeNumScheme; } + /*! + * \brief Get the kind of integration scheme (explicit or implicit) + * for the poisson/pressure correction equations. + * \note This value is obtained from the config file, and it is constant + * during the computation. + * \return Kind of integration scheme for the poisson/pressure correction equations. + */ + unsigned short GetKind_TimeIntScheme_Poisson(void) const { return Kind_TimeIntScheme_Poisson; } + /*! * \brief Get the kind of convective numerical scheme. * \note This is the information that the code will use, the method will diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 73266c7dacd..055f0a02c84 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -261,6 +261,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. */ }; @@ -336,6 +337,33 @@ static const MapType MatComp_Map = { MakePair("NEARLY_INCOMPRESSIBLE", STRUCT_COMPRESS::NEARLY_INCOMP) }; +/*! + * \brief Type of incompressible solver + */ +enum class ENUM_INCOMP_SYSTEM { + DENSITY_BASED, /*!< \brief Density-based. */ + PRESSURE_BASED, /*!< \brief Pressure-based. */ +}; +static const MapType Incomp_Map = { + MakePair("DENSITY_BASED", ENUM_INCOMP_SYSTEM::DENSITY_BASED) + MakePair("PRESSURE_BASED", ENUM_INCOMP_SYSTEM::PRESSURE_BASED) +}; + +/*! + * \brief Type of iteration + */ +enum class ENUM_PBITER { + SIMPLE, /*!< \brief SIMPLE algorithm. */ + SIMPLEC, /*!< \brief SIMPLE algorithm. */ + PISO, /*!< \brief PISO algorithm. */ +}; + +static const MapType PBIter_Map = { + MakePair("SIMPLE", ENUM_PBITER::SIMPLE) + MakePair("SIMPLEC", ENUM_PBITER::SIMPLEC) + MakePair("PISO", ENUM_PBITER::PISO) +}; + /*! * \brief Types of interpolators */ @@ -476,6 +504,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 { @@ -494,6 +523,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. */ }; @@ -825,7 +855,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) @@ -834,6 +865,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) }; @@ -860,7 +892,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) @@ -882,6 +915,7 @@ static const MapType Upwind_Map = { MakePair("SLAU2", UPWIND::SLAU2) MakePair("FDS", UPWIND::FDS) MakePair("LAX-FRIEDRICH", UPWIND::LAX_FRIEDRICH) + MakePair("UDS", UPWIND::UDS) }; /*! @@ -2694,6 +2728,7 @@ enum PERIODIC_QUANTITIES { PERIODIC_LIM_PRIM_1 , /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */ PERIODIC_LIM_PRIM_2 , /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */ PERIODIC_IMPLICIT , /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */ + PERIODIC_PRESSURE , /*!< \brief Corrected pressure communication. */ }; /*! @@ -2727,6 +2762,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. */ + PRESSURE_VAR , /*!< \brief Primitive variable (pressure) communication. */ + MASS_FUX , /*!< \brief Primitive variable (mass flux) communication. */ + MOM_COEFF , }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 0f7a5c6821f..ed2dd388b23 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1137,6 +1137,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, ENUM_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, ENUM_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*/ @@ -1880,6 +1885,8 @@ void CConfig::SetConfig_Options() { addEnumOption("TIME_DISCRE_HEAT", Kind_TimeIntScheme_Heat, Time_Int_Map, EULER_IMPLICIT); /* DESCRIPTION: Time discretization */ addEnumOption("TIMESTEP_HEAT", Kind_TimeStep_Heat, Heat_TimeStep_Map, MINIMUM); + /* DESCRIPTION: Time discretization */ + addEnumOption("TIME_DISCRE_POISSON", Kind_TimeIntScheme_Poisson, Time_Int_Map, EULER_IMPLICIT); /*!\par CONFIG_CATEGORY: Linear solver definition \ingroup Config*/ /*--- Options related to the linear solvers ---*/ @@ -1894,6 +1901,8 @@ 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); + /* DESCRIPTION: Maximum number of iterations of the poisson linear solver for the implicit formulation */ + addUnsignedLongOption("POISSON_LINEAR_SOLVER_ITER", Poisson_Linear_Solver_Iter, 10); /* DESCRIPTION: Fill in level for the ILU preconditioner */ addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0); /* DESCRIPTION: Use level scheduling for OMP parallelization of the ILU preconditioner */ @@ -1908,6 +1917,10 @@ 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 based flow corrections */ + addDoubleOption("RELAXATION_FACTOR_PBFLOW", Relaxation_Factor_PBFlow, 0.5); + /* DESCRIPTION: Relaxation of the Rhie Chow interpolation contribution in pressure based flow. */ + addDoubleOption("RELAXATION_FACTOR_RHIECHOW", RCFactor, 0.0); /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */ addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0); /* DESCRIPTION: Relaxation of the CHT coupling */ @@ -8805,6 +8818,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; @@ -8943,6 +8957,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(Kind_TimeIntScheme_Poisson); + } + break; + case MAIN_SOLVER::FEM_ELASTICITY: case MAIN_SOLVER::DISC_ADJ_FEM: if (val_system == RUNTIME_FEA_SYS) { diff --git a/TestCases/incomp_pressure_based/lid_driven_cavity/lid_driven_cavity.cfg b/TestCases/incomp_pressure_based/lid_driven_cavity/lid_driven_cavity.cfg new file mode 100644 index 00000000000..8774c960471 --- /dev/null +++ b/TestCases/incomp_pressure_based/lid_driven_cavity/lid_driven_cavity.cfg @@ -0,0 +1,272 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Incompressible and laminar lid driven cavity % +% Author: Akshay Koodly % +% Institution: University of Twente % +% File Version 7.0.7 "Blackbird % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +% Physical governing equations (EULER, NAVIER_STOKES, +% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, +% POISSON_EQUATION) + +SOLVER= INC_NAVIER_STOKES +% +% Incompressible solver (DENSITY_BASED, PRESSURE_BASED) +% DENSITY_BASED is chosen by default +KIND_INCOMP_SYSTEM= PRESSURE_BASED +% +% Gravity force (NO, YES) +GRAVITY_FORCE= NO +% +% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT) +MATH_PROBLEM= DIRECT +% +% Restart solution (NO, YES) +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +% Density model within the incompressible flow solver. +% Options are CONSTANT (default), BOUSSINESQ, or VARIABLE. If VARIABLE, +% an appropriate fluid model must be selected. +INC_DENSITY_MODEL= CONSTANT +% +% Solve the energy equation in the incompressible flow solver +INC_ENERGY_EQUATION = NO +% +% Initial density for incompressible flows (1.2886 kg/m^3 by default) +INC_DENSITY_INIT= 1.0 +% +% Reference density for incompressible flows (1.0 kg/m^3 by default) +INC_DENSITY_REF= 1.0 +% +% Reference velocity for incompressible flows (1.0 m/s by default) +INC_VELOCITY_REF= 1.0 +% +% Reference temperature for incompressible flows that include the +% energy equation (1.0 K by default) +INC_TEMPERATURE_REF = 1.0 +% Initial velocity for incompressible flows (1.0,0,0 m/s by default) +INC_VELOCITY_INIT= ( 1.0, 0.0, 0.0 ) + +% Non-dimensionalization scheme for incompressible flows. Options are +% INITIAL_VALUES (default), REFERENCE_VALUES, or DIMENSIONAL. +% INC_*_REF values are ignored unless REFERENCE_VALUES is chosen. +%INC_NONDIM= REFERENCE_VALUES +% +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY). +VISCOSITY_MODEL= CONSTANT_VISCOSITY +% +% Molecular Viscosity thsat would be constant (1.716E-5 by default) +MU_CONSTANT= 2.5E-3 +% +% +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +% Reference origin for moment computation +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +% +% Reference length for pitching, rolling, and yawing non-dimensional moment +REF_LENGTH= 1.0 +% +% Reference area for force coefficients (0 implies automatic calculation) +REF_AREA= 1.0 + +% ----------------------- BOUNDARY CONDITION DEFINITION -----------------------% +% +% Euler wall boundary marker(s) (NONE = no marker) +MARKER_HEATFLUX= ( Bottom,0.0 , Left, 0.0, Right, 0.0, Top, 0.0) +% +% Incompressible: (inlet marker, NULL, velocity magnitude, flow_direction_x, +% flow_direction_y, flow_direction_z, ... ) where flow_direction is +% a unit vector. +%MARKER_INLET= (Left, 0.0, 1.0, 1.0, 0.0, 0.0 ) +% +%INC_INLET_TYPE= VELOCITY_INLET +% +% Outlet boundary marker(s) (NONE = no marker) +%MARKER_OUTLET= ( Right, 0.0,Top, 0.0, Bottom, 0.0 ) +% +%INC_OUTLET_TYPE= PRESSURE_OUTLET, OPEN, OPEN +% Outlet boundary marker(s) (NONE = no marker) +%MARKER_OUTLET= ( Right, 0.0, Top, 0.0 ) +% +%INC_OUTLET_TYPE= PRESSURE_OUTLET,PRESSURE_OUTLET +% +% Symmetry boundary marker(s) (NONE = no marker) +%MARKER_EULER= ( Bottom) +% +%MARKER_SYM= (Bottom) +% Marker(s) of the surface that is going to be analyzed in detail (massflow, average pressure, distortion, etc) +%MARKER_ANALYZE = ( Left, Right, Top, Bottom ) + +% Method to compute the average value in MARKER_ANALYZE (AREA, MASSFLUX). +%MARKER_ANALYZE_AVERAGE = AREA + +%MARKER_PLOTTING= ( Bottom, Top ) +% +% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated +MARKER_MONITORING= ( Bottom, Top) + +% ----------------------- 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 +% +% +% Coordinates of the motion origin +SURFACE_MOTION_ORIGIN = 0.0 0.0 0.0 +% +% Angular velocity vector (rad/s) about the motion origin +SURFACE_ROTATION_RATE = 0.0 0.0 0.0 +SURFACE_TRANSLATION_RATE = 1.0 0.0 0.0 + +% ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------% +% +% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) +%NUM_METHOD_GRAD= GREEN_GAUSS +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +% +% Courant-Friedrichs-Lewy condition of the finest grid +CFL_NUMBER= 1.0 +% +% Adaptive CFL number (NO, YES) +CFL_ADAPT= NO +% +% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, +% CFL max value ) +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +% +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +% +% Number of total iterations +ITER= 10000 + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +% Multi-Grid Levels (0 = no multi-grid) +MGLEVEL= 0 +% +% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE) +MGCYCLE= V_CYCLE +% +% Multi-grid pre-smoothing level +MG_PRE_SMOOTH= ( 1, 2, 3, 3 ) +% +% Multi-grid post-smoothing level +MG_POST_SMOOTH= ( 1, 1, 1, 1 ) +% +% Jacobi implicit smoothing of the correction +MG_CORRECTION_SMOOTH= ( 1, 1, 1, 1 ) +% +% Damping factor for the residual restriction +MG_DAMP_RESTRICTION= 0.85 +% +% Damping factor for the correction prolongation +MG_DAMP_PROLONGATION= 0.85 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +% Linear solver for implicit formulations (BCGSTAB, FGMRES) +LINEAR_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS) +LINEAR_SOLVER_PREC= LU_SGS +% +% Minimum error of the linear solver for implicit formulations +LINEAR_SOLVER_ERROR= 1E-6 +% +% Max number of iterations of the linear solver for the implicit formulation +LINEAR_SOLVER_ITER= 20 +% +% Number of iterations for the poisson solver +POISSON_LINEAR_SOLVER_ITER= 50 +% +% Number of iterations for the flow solver +%LINEAR_SOLVER_ITER_FLOW= 25 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +% Convective numerical method (JST, JST_KE, JST_MAT, LAX-FRIEDRICH, CUSP, ROE, AUSM, +% AUSMPLUSUP, AUSMPLUSUP2, AUSMPWPLUS, HLLC, TURKEL_PREC, +% SW, MSW, FDS, SLAU, SLAU2, L2ROE, LMROE, UDS) +CONV_NUM_METHOD_FLOW= UDS +% +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_FLOW= YES +% +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_FLOW= NONE +% +% Coefficient for the limiter (smooth regions) +VENKAT_LIMITER_COEFF= 0.03 +% +% 2nd and 4th order artificial dissipation coefficients +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +% +% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +% Convective numerical method (SCALAR_UPWIND) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +% +% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. +% Required for 2nd order upwind schemes (NO, YES) +MUSCL_TURB= NO +% +% Slope limiter (VENKATAKRISHNAN, MINMOD) +SLOPE_LIMITER_TURB= VENKATAKRISHNAN +% +% Time discretization (EULER_IMPLICIT) +TIME_DISCRE_TURB= EULER_IMPLICIT +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% Convergence criteria (CAUCHY, RESIDUAL) +% +% CONV_CRITERIA= RESIDUAL +% +% Residual reduction (order of magnitude with respect to the initial value) +% RESIDUAL_REDUCTION= 7 +% +% Min value of the residual (log10 of the residual) +CONV_RESIDUAL_MINVAL= -10 +% +% Start Cauchy criteria at iteration number +CONV_STARTITER= 10 +% +% Number of elements to apply the criteria +CONV_CAUCHY_ELEMS= 50 +% +% Epsilon to control the series convergence +CONV_CAUCHY_EPS= 1E-6 +% +% Function to apply the criteria (LIFT, DRAG, SENS_GEOMETRY, SENS_MACH, +% DELTA_LIFT, DELTA_DRAG) +% CAUCHY_FUNC_FLOW= DRAG + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= square_65x65.su2 +MESH_FORMAT= SU2 +SOLUTION_FILENAME= solution_flow +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow +VOLUME_FILENAME= flow +SURFACE_FILENAME= surface_flow +OUTPUT_WRT_FREQ= 50 From 1a7e82a1d4cd1c45bd7ecc3aac4fc54f1a97366c Mon Sep 17 00:00:00 2001 From: Thijs Aalbers Date: Sat, 9 May 2026 22:30:56 +0200 Subject: [PATCH 2/4] Numerics and driver additions, WIP --- Common/include/option_structure.hpp | 8 +- SU2_CFD/include/numerics/CNumerics.hpp | 103 +++++ SU2_CFD/include/numerics/pbflow.hpp | 232 ++++++++++ SU2_CFD/include/numerics/poisson.hpp | 219 ++++++++++ SU2_CFD/src/drivers/CDriver.cpp | 139 ++++-- SU2_CFD/src/meson.build | 2 + SU2_CFD/src/numerics/CNumerics.cpp | 195 ++++++++- SU2_CFD/src/numerics/pbflow.cpp | 573 +++++++++++++++++++++++++ SU2_CFD/src/numerics/poisson.cpp | 421 ++++++++++++++++++ 9 files changed, 1853 insertions(+), 39 deletions(-) create mode 100644 SU2_CFD/include/numerics/pbflow.hpp create mode 100644 SU2_CFD/include/numerics/poisson.hpp create mode 100644 SU2_CFD/src/numerics/pbflow.cpp create mode 100644 SU2_CFD/src/numerics/poisson.cpp diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 055f0a02c84..c4f70e37eef 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. */ @@ -341,7 +341,7 @@ static const MapType MatComp_Map = { * \brief Type of incompressible solver */ enum class ENUM_INCOMP_SYSTEM { - DENSITY_BASED, /*!< \brief Density-based. */ + DENSITY_BASED, /*!< \brief Density-based. */ PRESSURE_BASED, /*!< \brief Pressure-based. */ }; static const MapType Incomp_Map = { @@ -353,8 +353,8 @@ static const MapType Incomp_Map = { * \brief Type of iteration */ enum class ENUM_PBITER { - SIMPLE, /*!< \brief SIMPLE algorithm. */ - SIMPLEC, /*!< \brief SIMPLE algorithm. */ + SIMPLE, /*!< \brief SIMPLE algorithm. */ + SIMPLEC, /*!< \brief SIMPLE algorithm. */ PISO, /*!< \brief PISO algorithm. */ }; diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 535e40cba96..9547164d92a 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -57,6 +57,7 @@ class CNumerics { su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl's number. */ su2double MassFlux; /*!< \brief Mass flux across edge. */ su2double + **Flux_Tensor, /*!< \brief Flux tensor (used for viscous and inviscid purposes). */ *Proj_Flux_Tensor; /*!< \brief Flux tensor projected in a direction. */ su2double **tau; /*!< \brief Viscous stress tensor. */ const su2double delta [3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; /*!< \brief Identity matrix. */ @@ -156,6 +157,12 @@ class CNumerics { TurbPsi_Grad_j, /*!< \brief Gradient of adjoint turbulent variables at point j. */ AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ + su2double + Poisson_Coeff_i, /*!< \brief Poisson coefficient at point i. */ + Poisson_Coeff_j; /*!< \brief Poisson coefficient at point j. */ + su2double + Poissonval_i, /*!< \brief Poisson value at point i. */ + Poissonval_j; /*!< \brief Poisson value at point j. */ su2double LocalGridLength_i; /*!< \brief Local grid length at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ @@ -180,6 +187,7 @@ class CNumerics { su2double StrainMag_i, StrainMag_j; /*!< \brief Strain rate magnitude. */ su2double Dissipation_i, Dissipation_j; /*!< \brief Dissipation. */ su2double Dissipation_ij; + su2double Source_Term; /*!< \brief Source term for poisson equation. */ su2double roughness_i = 0.0, /*!< \brief Roughness of the wall nearest to point i. */ roughness_j = 0.0; /*!< \brief Roughness of the wall nearest to point j. */ @@ -199,6 +207,7 @@ class CNumerics { su2double uq_delta_b; /*!< \brief Magnitude of perturbation */ su2double uq_urlx; /*!< \brief Under-relaxation factor for numerical stability */ bool uq_permute; /*!< \brief Flag for eigenvector permutation */ + su2double Face_Flux; /*!< \brief Face flux */ bool nemo; /*!< \brief Flag for NEMO problems */ @@ -1104,6 +1113,40 @@ class CNumerics { */ inline su2double GetDissipation() const { return Dissipation_ij; } + /*! + * \brief Set the Poisson coefficient + * \param[in] val_Poisson_i - Value of the Poisson coefficient at point i. + * \param[in] val_Poisson_j - Value of the Poisson coefficient at point j. + */ + inline void SetPoisson_Coeff(su2double val_Poisson_Coeff_i,su2double val_Poisson_Coeff_j){ + Poisson_Coeff_i = val_Poisson_Coeff_i; + Poisson_Coeff_j = val_Poisson_Coeff_j; + } + + /*! + * \brief Set the Poisson value + * \param[in] val_Poisson_i - Value of the Poisson variable at point i. + * \param[in] val_Poisson_j - Value of the Poisson variable at point j. + */ + inline void SetPoissonval(su2double val_Poisson_i,su2double val_Poisson_j){ + Poissonval_i = val_Poisson_i; + Poissonval_j = val_Poisson_j; + } + + /*! + * \brief Set the Poisson value + * \param[in] val_Poisson_i - Value of the Poisson variable at point i. + * \param[in] val_Poisson_j - Value of the Poisson variable at point j. + */ + inline void SetSourcePoisson(su2double val_Source_Term) { Source_Term = val_Source_Term; } + + /*! + * \brief Set the value of the momentum equation coefficients for Poisson eq. + * \param[in] val_Mom_Coeff_i - Value of the cross coefficient at point i. + * \param[in] val_Mom_Coeff_j - Value of the cross coefficient at point j. + */ + inline virtual void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { } + /*! * \brief Compute the projected inviscid flux vector. * \param[in] val_density - Pointer to the density. @@ -1177,6 +1220,21 @@ class CNumerics { su2double val_scale, su2double **val_Proj_Jac_Tensor) const; + /*! + * \brief Compute the projected inviscid flux vector for (pressure based) incompresible simulations + * \param[in] val_density - Pointer to the density. + * \param[in] val_velocity - Pointer to the velocity. + * \param[in] val_pressure - Pointer to the pressure. + * \param[in] val_betainc2 - Value of the artificial compresibility factor. + * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. + * \param[out] val_Proj_Flux - Pointer to the projected flux. + */ + void GetInviscidPBProjFlux(const su2double *val_density, + const su2double *val_velocity, + const su2double *val_pressure, + const su2double *val_normal, + su2double *val_Proj_Flux); + /*! * \brief Compute the low speed preconditioning matrix. * \param[in] val_density - Value of the density. @@ -1230,6 +1288,19 @@ class CNumerics { const su2double *val_Mean_SecVar, su2double **val_Jac_PC) const; + /*! + * \brief Compute the projection of the inviscid Jacobian matrices (Pressure based method). + * \param[in] val_density - Value of the density. + * \param[in] val_velocity - Pointer to the velocity. + * \param[in] val_cp - Value of the specific heat at constant pressure. + * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. + * \param[in] val_scale - Scale of the projection. + * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. + */ + void GetInviscidPBProjJac(const su2double val_density, const su2double *val_velocity, + const su2double *val_normal, const su2double val_scale, + su2double **val_Proj_Jac_tensor); + /*! * \overload * \brief Computation of the matrix P for a generic fluid model @@ -1854,6 +1925,38 @@ class CNumerics { * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport */ inline bool GetBoundedScalar() const { return bounded_scalar;} + + /*! + * \brief Compute the projection of the viscous Jacobian matrices. + * \param[in] val_laminar_viscosity - Value of the laminar viscosity. + * \param[in] val_eddy_viscosity - Value of the eddy viscosity. + * \param[in] val_dist_ij - Distance between the points. + * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. + * \param[in] val_dS - Area of the face between two nodes. + * \param[out] val_Proj_Jac_Tensor_i - Pointer to the projected viscous Jacobian at point i. + * \param[out] val_Proj_Jac_Tensor_j - Pointer to the projected viscous Jacobian at point j. + */ + void GetViscousPBIncProjJacs(const su2double val_laminar_viscosity, + const su2double val_eddy_viscosity, + const su2double val_dist_ij, + const su2double *val_normal, + const su2double val_dS, + su2double **val_Proj_Jac_Tensor_i, + su2double **val_Proj_Jac_Tensor_j); + /*! + * \brief Compute the projection of the viscous fluxes into a direction (pressure based method). + * \param[in] val_primvar - Primitive variables. + * \param[in] val_gradprimvar - Gradient of the primitive variables. + * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. + * \param[in] val_laminar_viscosity - Laminar viscosity. + * \param[in] val_eddy_viscosity - Eddy viscosity. + */ + void GetViscousPBIncProjFlux(const su2double *val_primvar, + su2double **val_gradprimvar, + const su2double *val_normal, + const su2double val_laminar_viscosity, + const su2double val_eddy_viscosity, + const su2double val_turb_ke); }; /*! diff --git a/SU2_CFD/include/numerics/pbflow.hpp b/SU2_CFD/include/numerics/pbflow.hpp new file mode 100644 index 00000000000..637db426d1f --- /dev/null +++ b/SU2_CFD/include/numerics/pbflow.hpp @@ -0,0 +1,232 @@ +/*! + * \file pbflow.hpp + * \brief Delarations of numerics classes for heat transfer problems. + * \author F. Palacios, T. Economon + * \version 8.0.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2020, 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 CUpwPB_Flow + * \brief Class for solving an approximate Riemann solver of Roe for the pressure based incompressible flow equations. + * \ingroup ConvDiscr + * \author F. Palacios + */ +class CUpwPB_Flow : public CNumerics { +private: + bool implicit, dynamic_grid; + bool gravity; + su2double Froude, Upw_i, Upw_j; + su2double *Diff_U; + su2double *Velocity_i, *Velocity_j, *MeanVelocity, *Velocity_upw; + su2double *ProjFlux_i, *ProjFlux_j; + su2double Proj_ModJac_Tensor_ij, Pressure_i, + Pressure_j, MeanDensity, MeanSoundSpeed, MeanPressure, MeanBetaInc2, + ProjVelocity, FaceVel, Face_Flux; + unsigned short iDim, iVar, jVar, kVar; + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + su2double **Jacobian_upw = nullptr; + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CUpwPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CUpwPB_Flow(void); + + /*! + * \brief Compute the Roe's flux between two nodes i and j. + * \param[out] val_residual - Pointer to the total 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) override; + + /*! + * \brief Set the value of face velocity. This is used as a proxy for massflux at a face. + * \param[in] val_FaceVel. + */ + inline void SetFaceVel(su2double val_FaceVel){ FaceVel = val_FaceVel; } +}; + + +/*! + * \class CCentLaxArtComp_Flow + * \brief Class for computing the Lax-Friedrich centered scheme (artificial compressibility). + * \ingroup ConvDiscr + * \author F. Palacios + */ +class CCentPB_Flow : public CNumerics { +private: + bool implicit, dynamic_grid; + bool gravity; + su2double Froude, Upw_i, Upw_j; + su2double *Diff_U; + su2double *Velocity_i, *Velocity_j, *MeanVelocity, *Velocity_upw; + su2double *ProjFlux_i, *ProjFlux_j; + su2double *Lambda, *Epsilon; + su2double **P_Tensor, **invP_Tensor,**val_Jacobian_upw; + su2double Proj_ModJac_Tensor_ij, Pressure_i, + Pressure_j, MeanDensity, MeanSoundSpeed, MeanPressure, MeanBetaInc2, + ProjVelocity, FaceVel; + unsigned short iDim, iVar, jVar, kVar; + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + su2double **Jacobian_upw = nullptr; + +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. + */ + CCentPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CCentPB_Flow(void); + + /*! + * \brief Compute the flow residual using a Lax method. + * \param[out] val_resconv - Pointer to the convective residual. + * \param[out] val_resvisc - Pointer to the artificial viscosity 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) override; +}; + +/*! + * \class CAvgGradInc_Flow + * \brief Class for computing viscous term using an average of gradients. + * \ingroup ViscDiscr + * \author A. Bueno, F. Palacios, T. Economon + */ +class CAvgGradPBInc_Flow : public CNumerics { +private: + unsigned short iDim, iVar, jVar; /*!< \brief Iterators in dimension an variable. */ + su2double *Mean_PrimVar, /*!< \brief Mean primitive variables. */ + *PrimVar_i, *PrimVar_j; /*!< \brief Primitives variables at point i and j. */ + su2double **Mean_GradPrimVar, /*!< \brief Mean value of the gradient. */ + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, /*!< \brief Mean value of the viscosity. */ + Mean_turb_ke, /*!< \brief Mean value of the turbulent kinetic energy. */ + dist_ij, /*!< \brief Length of the edge and face. */ + proj_vector_ij; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ + bool implicit; /*!< \brief Implicit calculus. */ + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + +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. + */ + CAvgGradPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CAvgGradPBInc_Flow(void); + /*! + * \brief Compute the viscous flow residual using an average of gradients. + * \param[out] val_residual - Pointer to the total 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) override; +}; + + +/*! + * \class CAvgGradCorrectedPBInc_Flow + * \brief Class for computing viscous term using an average of gradients with correction (incompressible). + * \ingroup ViscDiscr + */ +class CAvgGradCorrectedPBInc_Flow : public CNumerics { +private: + unsigned short iDim, iVar, jVar; /*!< \brief Iterators in dimension an variable. */ + su2double *Mean_PrimVar; /*!< \brief Mean primitive variables. */ + su2double *PrimVar_i, *PrimVar_j, /*!< \brief Primitives variables at point i and 1. */ + *Edge_Vector, /*!< \brief Vector form point i to point j. */ + **Mean_GradPrimVar, *Proj_Mean_GradPrimVar_Edge, /*!< \brief Mean value of the gradient. */ + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, /*!< \brief Mean value of the viscosity. */ + Mean_turb_ke, /*!< \brief Mean value of the turbulent kinetic energy. */ + dist_ij_2, /*!< \brief Length of the edge and face. */ + proj_vector_ij; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ + bool implicit; /*!< \brief Implicit calculus. */ + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + +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. + */ + CAvgGradCorrectedPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CAvgGradCorrectedPBInc_Flow(void); + + /*! + * \brief Compute the viscous flow residual using an average of gradients with correction. + * \param[out] val_residual - Pointer to the total 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) override; +}; diff --git a/SU2_CFD/include/numerics/poisson.hpp b/SU2_CFD/include/numerics/poisson.hpp new file mode 100644 index 00000000000..21fe28d60b0 --- /dev/null +++ b/SU2_CFD/include/numerics/poisson.hpp @@ -0,0 +1,219 @@ +/*! + * \file heat.hpp + * \brief Delarations of numerics classes for heat transfer problems. + * \author F. Palacios, T. Economon + * \version 8.0.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2020, 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 CAvgGrad_Poisson + * \brief Class for computing viscous term using average of gradients without correction (Poisson equation). + * \ingroup ViscDiscr + */ +class CAvgGrad_Poisson : public CNumerics { +private: + + su2double *Edge_Vector; + bool implicit,direct; + su2double **Mean_GradPoissonVar; + su2double *Mom_Coeff_i,*Mom_Coeff_j; + su2double *Proj_Mean_GradPoissonVar_Normal, *Proj_Mean_GradPoissonVar_Corrected; + su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean ; + unsigned short iVar, iDim; + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CAvgGrad_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CAvgGrad_Poisson(void); + + /*! + * \brief Compute the viscous heat residual using an average of gradients with correction. + * \param[out] val_residual - Pointer to the total residual. + * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[out] 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) override; + + + inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; + Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; + } + } +}; + +/*! + * \class CAvgGradCorrected_Poisson + * \brief Class for computing viscous term using average of gradients with correction (Poisson equation). + * \ingroup ViscDiscr + */ + +class CAvgGradCorrected_Poisson : public CNumerics { +private: + + su2double *Edge_Vector; + bool implicit,direct; + su2double **Mean_GradPoissonVar; + su2double *Mom_Coeff_i,*Mom_Coeff_j; + su2double *Proj_Mean_GradPoissonVar_Kappa, *Proj_Mean_GradPoissonVar_Edge, *Proj_Mean_GradPoissonVar_Corrected; + su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean; + unsigned short iVar, iDim; + + su2double *Flux = nullptr; + su2double **Jacobian_i = nullptr; + su2double **Jacobian_j = nullptr; + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CAvgGradCorrected_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CAvgGradCorrected_Poisson(void); + + /*! + * \brief Compute the viscous residual using an average of gradients with correction. + * \param[out] val_residual - Pointer to the total residual. + * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[out] 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) override; + + + inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; + Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; + } + } +}; + +/*! + * \class CAvgGrad_Poisson + * \brief Class for computing viscous term using average of gradients without correction (Poisson equation). + * \ingroup ViscDiscr + */ +class CPressure_Poisson : public CNumerics { +private: + + su2double *Edge_Vector; + bool implicit,direct; + su2double **Mean_GradPoissonVar; + su2double *Proj_Mean_GradPoissonVar_Kappa, *Proj_Mean_GradPoissonVar_Edge, *Proj_Mean_GradPoissonVar_Corrected; + su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean; + su2double *Mom_Coeff_i,*Mom_Coeff_j; + unsigned short iVar, iDim; + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CPressure_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CPressure_Poisson(void); + + /*! + * \brief Compute the viscous heat residual using an average of gradients with correction. + * \param[out] val_residual - Pointer to the total residual. + * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[out] Jacobian_j - Jacobian of the numerical method at node j (implicit computation). + * \param[in] config - Definition of the particular problem. + */ + void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config); + + + inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; + Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; + } + } +}; + +/*! + * \class CSource_Poisson + * \brief Class for source term of the Poisson equation. + * \ingroup SourceDiscr + */ +class CSource_PoissonFVM : public CNumerics { +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Name of the input config file + * + */ + CSource_PoissonFVM(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + + /*! + * \brief Residual for source term integration. + * \param[out] val_residual - Pointer to the total residual. + * \param[out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[in] config - Definition of the particular problem. + */ + void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CSource_PoissonFVM(void); +}; + diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 5bdeb9ce1af..c95474c9ef7 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -56,6 +56,8 @@ #include "../../include/numerics/template.hpp" #include "../../include/numerics/radiation.hpp" #include "../../include/numerics/heat.hpp" +#include "../../include/numerics/pbflow.hpp" +#include "../../include/numerics/poisson.hpp" #include "../../include/numerics/flow/convection/roe.hpp" #include "../../include/numerics/flow/convection/fds.hpp" #include "../../include/numerics/flow/convection/fvs.hpp" @@ -1447,13 +1449,16 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver nVar_Adj_Flow = 0, nVar_Adj_Turb = 0, nVar_FEM = 0, - nVar_Rad = 0; + nVar_Rad = 0, + nVar_Poisson = 0; numerics = new CNumerics***[config->GetnMGLevels()+1] (); 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() == ENUM_INCOMP_SYSTEM::PRESSURE_BASED); + bool poisson = (config->GetKind_Incomp_System() == ENUM_INCOMP_SYSTEM::PRESSURE_BASED); bool roe_low_dissipation = (config->GetKind_RoeLowDiss() != NO_ROELOWDISS); /*--- Initialize some useful booleans ---*/ @@ -1567,6 +1572,8 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver if (fem) nVar_FEM = solver[MESH_0][FEA_SOL]->GetnVar(); + if (poisson) nVar_Poisson = solver[MESH_0][POISSON_SOL]->GetnVar(); + if (config->AddRadiation()) nVar_Rad = solver[MESH_0][RAD_SOL]->GetnVar(); /*--- Number of variables for adjoint problem ---*/ @@ -1654,22 +1661,37 @@ 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 CCentPB_Flow(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 CCentPB_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 CUpwPB_Flow(nDim, nVar_Flow, config); + + } } break; case SPACE_UPWIND : @@ -1776,17 +1798,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 CUpwPB_Flow(nDim, nVar_Flow, config); + numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwPB_Flow(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; @@ -1823,14 +1860,25 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } } if (incompressible) { - /*--- Incompressible flow, use preconditioning method ---*/ - numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, true, config); - for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + if (!pressure_based) { + /*--- Incompressible flow, use preconditioning method ---*/ + numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, true, config); + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); - /*--- Definition of the boundary condition method ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + /*--- Definition of the boundary condition method ---*/ + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); + } else { + /*--- Incompressible flow, use pressure based method ---*/ + numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradCorrectedPBInc_Flow(nDim, nVar_Flow, config); + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); + + /*--- Definition of the boundary condition method ---*/ + for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); + } } /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ @@ -2035,8 +2083,13 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver if (turbulent) { if (incompressible) - InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + if (!pressure_based) + InstantiateTurbulentNumerics >(nVar_Turb, offset, config, solver[MESH_0][TURB_SOL], numerics); + // else + // InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + // solver[MESH_0][TURB_SOL], numerics); + else if (NEMO_ns) InstantiateTurbulentNumerics >(nVar_Turb, offset, config, solver[MESH_0][TURB_SOL], numerics); @@ -2094,6 +2147,24 @@ 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 CAvgGradCorrected_Poisson(nDim, nVar_Poisson, config); + + for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + numerics[iMGlevel][POISSON_SOL][visc_term] = new CAvgGrad_Poisson(nDim, nVar_Poisson, config); + + /*--- 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_Poisson(nDim, nVar_Poisson, config); + + /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ + numerics[iMGlevel][POISSON_SOL][source_first_term] = new CSource_PoissonFVM(nDim, nVar_Poisson, config); + numerics[iMGlevel][POISSON_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Poisson, config); + } + } + /*--- Solver definition for the radiation model problem ---*/ if (config->AddRadiation()) { diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3b40822a34f..6b754aec007 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -122,6 +122,8 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/template.cpp', 'numerics/radiation.cpp', + 'numerics/poisson.cpp', + 'numerics/pbflow.cpp', 'numerics/flow/convection/roe.cpp', 'numerics/flow/convection/fds.cpp', 'numerics/flow/convection/fvs.cpp', diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index bd8728e0e83..6ff3efcc9aa 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -44,7 +44,7 @@ CNumerics::CNumerics() { CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) { - unsigned short iDim; + unsigned short iDim, iVar; nDim = val_nDim; nVar = val_nVar; @@ -55,6 +55,10 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Gas_Constant = config->GetGas_ConstantND(); energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); + Flux_Tensor = new su2double* [nVar]; + for (iVar = 0; iVar < (nVar); iVar++) + Flux_Tensor[iVar] = new su2double [nDim] (); + tau = new su2double* [nDim]; for (iDim = 0; iDim < nDim; iDim++) tau[iDim] = new su2double [nDim] (); @@ -80,6 +84,12 @@ CNumerics::~CNumerics() { // visc delete [] Proj_Flux_Tensor; + if (Flux_Tensor) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) + delete [] Flux_Tensor[iVar]; + delete [] Flux_Tensor; + } + if (tau) { for (unsigned short iDim = 0; iDim < nDim; iDim++) delete [] tau[iDim]; @@ -1575,3 +1585,186 @@ su2double CNumerics::GetRoe_Dissipation(const su2double Dissipation_i, } return Dissipation_ij; } + +void CNumerics::GetInviscidPBProjFlux(const su2double *val_density, + const su2double *val_velocity, + const su2double *val_pressure, + const su2double *val_normal, + su2double *val_Proj_Flux) { + su2double rhou, rhov, rhow; + + if (nDim == 2) { + rhou = (*val_density)*val_velocity[0]; + rhov = (*val_density)*val_velocity[1]; + + val_Proj_Flux[0] = (rhou*val_velocity[0])*val_normal[0] + rhou*val_velocity[1]*val_normal[1]; + val_Proj_Flux[1] = rhov*val_velocity[0]*val_normal[0] + (rhov*val_velocity[1])*val_normal[1]; + } + else { + rhou = (*val_density)*val_velocity[0]; + rhov = (*val_density)*val_velocity[1]; + rhow = (*val_density)*val_velocity[2]; + + val_Proj_Flux[0] = (rhou*val_velocity[0])*val_normal[0] + rhou*val_velocity[1]*val_normal[1] + rhou*val_velocity[2]*val_normal[2]; + val_Proj_Flux[1] = rhov*val_velocity[0]*val_normal[0] + (rhov*val_velocity[1])*val_normal[1] + rhov*val_velocity[2]*val_normal[2]; + val_Proj_Flux[2] = rhow*val_velocity[0]*val_normal[0] + rhow*val_velocity[1]*val_normal[1] + (rhow*val_velocity[2])*val_normal[2]; + } +} + +void CNumerics::GetInviscidPBProjJac(const su2double val_density, const su2double *val_velocity, const su2double *val_normal, + const su2double val_scale, su2double **val_Proj_Jac_Tensor) { + + unsigned short iDim; + su2double proj_vel; + + proj_vel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + proj_vel += val_velocity[iDim]*val_normal[iDim]; + + if (nDim == 2) { + val_Proj_Jac_Tensor[0][0] = val_scale*val_density*(val_velocity[0]*val_normal[0] + proj_vel); + val_Proj_Jac_Tensor[0][1] = val_scale*val_density*val_velocity[0]*val_normal[1]; + + val_Proj_Jac_Tensor[1][0] = val_scale*val_density*val_velocity[1]*val_normal[0]; + val_Proj_Jac_Tensor[1][1] = val_scale*val_density*(val_velocity[1]*val_normal[1] + proj_vel); + } + else { + val_Proj_Jac_Tensor[0][0] = val_scale*val_density*(proj_vel+val_velocity[0]*val_normal[0]); + val_Proj_Jac_Tensor[0][1] = val_scale*val_density*(val_velocity[0]*val_normal[1]); + val_Proj_Jac_Tensor[0][2] = val_scale*val_density*(val_velocity[0]*val_normal[2]); + + val_Proj_Jac_Tensor[1][0] = val_scale*val_density*(val_velocity[1]*val_normal[0]); + val_Proj_Jac_Tensor[1][1] = val_scale*val_density*(proj_vel+val_velocity[1]*val_normal[1]); + val_Proj_Jac_Tensor[1][2] = val_scale*val_density*(val_velocity[1]*val_normal[2]); + + val_Proj_Jac_Tensor[2][0] = val_scale*val_density*(val_velocity[2]*val_normal[0]); + val_Proj_Jac_Tensor[2][1] = val_scale*val_density*(val_velocity[2]*val_normal[1]); + val_Proj_Jac_Tensor[2][2] = val_scale*val_density*(proj_vel+val_velocity[2]*val_normal[2]); + + } +} + +void CNumerics::GetViscousPBIncProjFlux(const su2double *val_primvar, + su2double **val_gradprimvar, + const su2double *val_normal, + const su2double val_laminar_viscosity, + const su2double val_eddy_viscosity, + const su2double val_turb_ke) { + + unsigned short iVar, iDim, jDim; + su2double total_viscosity, div_vel, Density; + + Density = val_primvar[nDim+1]; + + total_viscosity = (val_laminar_viscosity + val_eddy_viscosity); + + /*--- The full stress tensor is needed for variable density, as nabla.u != 0 ---*/ + /*--- Note: Gradients are computed as [iDim][jDim] and not [iDim+1][jDim] because + *--- the mean gradient passed here contains only velocities and no pressure. The + *--- mean gradient is computed as PrimVar_Grad[iVar+1][iDim], so no need to repeat.---*/ + + div_vel = 0.0; + for (iDim = 0 ; iDim < nDim; iDim++) + div_vel += val_gradprimvar[iDim][iDim]; + + for (iDim = 0 ; iDim < nDim; iDim++) + for (jDim = 0 ; jDim < nDim; jDim++) + tau[iDim][jDim] = (total_viscosity*(val_gradprimvar[jDim][iDim] + + val_gradprimvar[iDim][jDim] ) + -TWO3*total_viscosity*div_vel*delta[iDim][jDim] + -TWO3*Density*val_turb_ke*delta[iDim][jDim]); + + /*--- Gradient of primitive variables -> [Pressure vel_x vel_y vel_z Density] ---*/ + + if (nDim == 2) { + Flux_Tensor[0][0] = tau[0][0]; + Flux_Tensor[1][0] = tau[0][1]; + + Flux_Tensor[0][1] = tau[1][0]; + Flux_Tensor[1][1] = tau[1][1]; + + } else { + + Flux_Tensor[0][0] = tau[0][0]; + Flux_Tensor[1][0] = tau[0][1]; + Flux_Tensor[2][0] = tau[0][2]; + + Flux_Tensor[0][1] = tau[1][0]; + Flux_Tensor[1][1] = tau[1][1]; + Flux_Tensor[2][1] = tau[1][2]; + + Flux_Tensor[0][2] = tau[2][0]; + Flux_Tensor[1][2] = tau[2][1]; + Flux_Tensor[2][2] = tau[2][2]; + + } + + for (iVar = 0; iVar < nVar; iVar++) { + Proj_Flux_Tensor[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Proj_Flux_Tensor[iVar] += Flux_Tensor[iVar][iDim] * val_normal[iDim]; + } + +} + +void CNumerics::GetViscousPBIncProjJacs(const su2double val_laminar_viscosity, + const su2double val_eddy_viscosity, + const su2double val_dist_ij, + const su2double *val_normal, + const su2double val_dS, + su2double **val_Proj_Jac_Tensor_i, su2double **val_Proj_Jac_Tensor_j) { + unsigned short iDim, iVar, jVar; + + su2double theta = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + theta += val_normal[iDim]*val_normal[iDim]; + + su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; + su2double factor = total_viscosity/(val_dist_ij)*val_dS; + + if (nDim == 3) { + su2double thetax = theta + val_normal[0]*val_normal[0]/3.0; + su2double thetay = theta + val_normal[1]*val_normal[1]/3.0; + su2double thetaz = theta + val_normal[2]*val_normal[2]/3.0; + + su2double etax = val_normal[1]*val_normal[2]/3.0; + su2double etay = val_normal[0]*val_normal[2]/3.0; + su2double etaz = val_normal[0]*val_normal[1]/3.0; + + val_Proj_Jac_Tensor_i[0][0] = -factor*thetax; + val_Proj_Jac_Tensor_i[0][1] = -factor*etaz; + val_Proj_Jac_Tensor_i[0][2] = -factor*etay; + + val_Proj_Jac_Tensor_i[1][0] = -factor*etaz; + val_Proj_Jac_Tensor_i[1][1] = -factor*thetay; + val_Proj_Jac_Tensor_i[1][2] = -factor*etax; + + val_Proj_Jac_Tensor_i[2][0] = -factor*etay; + val_Proj_Jac_Tensor_i[2][1] = -factor*etax; + val_Proj_Jac_Tensor_i[2][2] = -factor*thetaz; + + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + val_Proj_Jac_Tensor_j[iVar][jVar] = -val_Proj_Jac_Tensor_i[iVar][jVar]; + + } + + if (nDim == 2) { + su2double thetax = theta + val_normal[0]*val_normal[0]/3.0; + su2double thetay = theta + val_normal[1]*val_normal[1]/3.0; + su2double etaz = val_normal[0]*val_normal[1]/3.0; + + + val_Proj_Jac_Tensor_i[0][0] = -factor*thetax; + val_Proj_Jac_Tensor_i[0][1] = -factor*etaz; + + val_Proj_Jac_Tensor_i[1][0] = -factor*etaz; + val_Proj_Jac_Tensor_i[1][1] = -factor*thetay; + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + val_Proj_Jac_Tensor_j[iVar][jVar] = -val_Proj_Jac_Tensor_i[iVar][jVar]; + } + +} diff --git a/SU2_CFD/src/numerics/pbflow.cpp b/SU2_CFD/src/numerics/pbflow.cpp new file mode 100644 index 00000000000..7a940bea026 --- /dev/null +++ b/SU2_CFD/src/numerics/pbflow.cpp @@ -0,0 +1,573 @@ +/*! + * \file pbflow.cpp + * \brief This file contains the numerical methods for incompressible flow. + * \author F. Palacios, T. Economon + * \version 6.0.1 "Falcon" + * + * The current SU2 release has been coordinated by the + * SU2 International Developers Society + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2018, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * 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/pbflow.hpp" + +CUpwPB_Flow::CUpwPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + gravity = config->GetGravityForce(); + Froude = config->GetFroude(); + dynamic_grid = config->GetDynamic_Grid(); + + Diff_U = new su2double [nVar]; + Velocity_i = new su2double [nDim]; + Velocity_j = new su2double [nDim]; + Velocity_upw = new su2double [nDim]; + MeanVelocity = new su2double [nDim]; + Flux = new su2double [nDim]; + ProjFlux_i = new su2double [nVar]; + ProjFlux_j = new su2double [nVar]; + Jacobian_i = new su2double* [nVar]; + Jacobian_j = new su2double* [nVar]; + Jacobian_upw = new su2double* [nVar]; + + + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nVar]; + Jacobian_j[iVar] = new su2double [nVar]; + Jacobian_upw[iVar] = new su2double [nVar]; + } + +} + +CUpwPB_Flow::~CUpwPB_Flow(void) { + + delete [] Diff_U; + delete [] Velocity_i; + delete [] Velocity_j; + delete [] Velocity_upw; + delete [] MeanVelocity; + delete [] Flux; + delete [] ProjFlux_i; + delete [] ProjFlux_j; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + delete [] Jacobian_upw[iVar]; + } + delete [] Jacobian_i; + delete [] Jacobian_j; + delete [] Jacobian_upw; + +} + +CNumerics::ResidualType<> CUpwPB_Flow::ComputeResidual(const CConfig *config) { + + + su2double MeanDensity, Flux0, Flux1, MeanPressure, Area, FF, Vel0, Vel1, ProjGridVelFlux = 0.0; + + + /*--- Primitive variables at point i and j ---*/ + Pressure_i = V_i[0]; Pressure_j = V_j[0]; + DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; + MeanDensity = 0.5*(DensityInc_i + DensityInc_j); + MeanPressure = 0.5*(Pressure_i + Pressure_j); + + Area = 0.0; + for(iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + + /*--- (rho*u_i) ---*/ + Face_Flux = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = V_i[iDim+1]; + Velocity_j[iDim] = V_j[iDim+1]; + MeanVelocity[iDim] = 0.5*(Velocity_i[iDim] + Velocity_j[iDim]); + Face_Flux += MeanDensity*MeanVelocity[iDim]*Normal[iDim]; + } + + if (dynamic_grid) { + ProjGridVelFlux = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVelFlux += 0.5*MeanDensity*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + } + Face_Flux -= ProjGridVelFlux; + } + + /*--- Find upwind direction. ---*/ + Flux0 = 0.5*(Face_Flux + fabs(Face_Flux)) ; + Flux1 = 0.5*(Face_Flux - fabs(Face_Flux)) ; + + Upw_i = round(fabs(Flux0/(fabs(Face_Flux)+EPS))); + Upw_j = round(fabs(Flux1/(fabs(Face_Flux)+EPS))); + + /*--- Find flux. ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + Flux[iVar] = Flux0*V_i[iVar+1] + Flux1*V_j[iVar+1]; + Velocity_upw[iVar] = Upw_i*V_i[iVar+1] + Upw_j*V_j[iVar+1]; + if (dynamic_grid) Velocity_upw[iVar] -= (Upw_i*GridVel_i[iVar] + Upw_j*GridVel_j[iVar]); + } + + if (implicit) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_j[iVar][jVar] = 0.0; + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_upw[iVar][jVar] = 0.0; + } + + GetInviscidPBProjJac(MeanDensity, Velocity_upw, Normal, 0.5, Jacobian_upw); + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = Upw_i*Jacobian_upw[iVar][jVar]; + Jacobian_j[iVar][jVar] = Upw_j*Jacobian_upw[iVar][jVar]; + } + } + + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); +} + +CCentPB_Flow::CCentPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + gravity = config->GetGravityForce(); + Froude = config->GetFroude(); + dynamic_grid = config->GetDynamic_Grid(); + + Diff_U = new su2double [nVar]; + Velocity_i = new su2double [nDim]; + Velocity_j = new su2double [nDim]; + Velocity_upw = new su2double [nDim]; + MeanVelocity = new su2double [nDim]; + Flux = new su2double [nDim]; + ProjFlux_i = new su2double [nVar]; + ProjFlux_j = new su2double [nVar]; + Jacobian_i = new su2double* [nVar]; + Jacobian_j = new su2double* [nVar]; + Jacobian_upw = new su2double* [nVar]; + + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nVar]; + Jacobian_j[iVar] = new su2double [nVar]; + Jacobian_upw[iVar] = new su2double [nVar]; + } +} + +CCentPB_Flow::~CCentPB_Flow(void) { + + delete [] Diff_U; + delete [] Velocity_i; + delete [] Velocity_j; + delete [] Velocity_upw; + delete [] MeanVelocity; + delete [] Flux; + delete [] ProjFlux_i; + delete [] ProjFlux_j; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + delete [] Jacobian_upw[iVar]; + } + delete [] Jacobian_i; + delete [] Jacobian_j; + delete [] Jacobian_upw; + +} + +CNumerics::ResidualType<> CCentPB_Flow::ComputeResidual(const CConfig *config) { + /*--- To do list + * 1. Find upwind direction, upwind node and downwind node + * 2. Compute normalised variable for the upwind node using the gradient of the upwind node + * 3. Find face velocity using central difference scheme + * 4. Use the ULTIMATE limiter to find the adjusted face velocity + * 5. Find residual as FaceFlux * AdjustedFaceVelocity + * */ + + + su2double MeanDensity, Flux0, Flux1, MeanPressure, Area, FF, Vel0, Vel1, ProjGridVelFlux = 0.0; + su2double dissipation, kappa=0.15, ProjVel_i = 0.0, ProjVel_j = 0.0; + + /*--- Conservative variables at point i and j ---*/ + + Pressure_i = V_i[0]; Pressure_j = V_j[0]; + DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; + + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = V_i[iDim+1]; + Velocity_j[iDim] = V_j[iDim+1]; + } + + Area = 0.0; + for(iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + + Face_Flux = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = V_i[iDim+1]; + Velocity_j[iDim] = V_j[iDim+1]; + MeanVelocity[iDim] = 0.5*(Velocity_i[iDim] + Velocity_j[iDim]); + Face_Flux += MeanDensity*MeanVelocity[iDim]*Normal[iDim]; + ProjVel_i += Velocity_i[iDim]*Normal[iDim]; + ProjVel_j += Velocity_j[iDim]*Normal[iDim]; + } + + if (dynamic_grid) { + ProjGridVelFlux = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVelFlux += 0.5*MeanDensity*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + } + Face_Flux -= ProjGridVelFlux; + } + + Flux0 = 0.5*(Face_Flux + fabs(Face_Flux)) ; + Flux1 = 0.5*(Face_Flux - fabs(Face_Flux)) ; + + Upw_i = round(fabs(Flux0/(fabs(Face_Flux)+EPS))); + Upw_j = round(fabs(Flux1/(fabs(Face_Flux)+EPS))); + + for (iVar = 0; iVar < nVar; iVar++) { + Velocity_upw[iVar] = Upw_i*V_i[iVar+1] + Upw_j*V_j[iVar+1]; + if (dynamic_grid) Velocity_upw[iVar] -= (Upw_i*GridVel_i[iVar] + Upw_j*GridVel_j[iVar]); + Flux[iVar] = Face_Flux*MeanVelocity[iVar]; + } + + //Dissipation + su2double lambda_i = 0.0, lambda_j = 0.0; + lambda_i = 2.0*abs(ProjVel_i); + lambda_j = 2.0*abs(ProjVel_j); + + su2double lambda_mean = 0.5*(lambda_i + lambda_j); + if (lambda_mean < EPS) { + lambda_i = 2.0*abs(config->GetVelocity_Ref()*Area); + lambda_j = 2.0*abs(config->GetVelocity_Ref()*Area); + lambda_mean = abs(config->GetVelocity_Ref()*Area); + } + + su2double Param_p = 0.3, SF=0.0; + su2double Phi_i = pow((lambda_i/(4.0*lambda_mean)),Param_p); + su2double Phi_j = pow((lambda_j/(4.0*lambda_mean)),Param_p); + + if ((Phi_i + Phi_j) != 0.0) + SF = 4.0*Phi_i*Phi_j/(Phi_i + Phi_j); + + su2double sc0 = 3.0*(Neighbor_i + Neighbor_j)/(Neighbor_i*Neighbor_j); + su2double E_0 = kappa*sc0*2.0/3.0; + + su2double diss[MAXNDIM]; + for (iDim = 0; iDim < nDim; iDim++) + diss[iDim] = E_0*(DensityInc_i*Velocity_i[iDim] - DensityInc_j*Velocity_j[iDim])*SF*lambda_mean; + + for (iVar = 0; iVar < nVar; iVar++) + Flux[iVar] += diss[iVar]; + + /*--- For implicit schemes, compute jacobians based on the upwind scheme for stability issues. (See ANSYS user guide) ---*/ + if (implicit) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_j[iVar][jVar] = 0.0; + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_upw[iVar][jVar] = 0.0; + } + + GetInviscidPBProjJac(MeanDensity, MeanVelocity, Normal, 0.5, Jacobian_upw); + //GetInviscidPBProjJac(&DensityInc_i, Velocity_upw, Normal, 0.5, val_Jacobian_upw); + /*GetInviscidPBProjJac(&DensityInc_i, Velocity_i, Normal, 0.5, Jacobian_i); + GetInviscidPBProjJac(&DensityInc_i, Velocity_j, Normal, 0.5, Jacobian_j);*/ + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = Upw_i*Jacobian_upw[iVar][jVar]; + Jacobian_j[iVar][jVar] = Upw_j*Jacobian_upw[iVar][jVar]; + } + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar][iVar] += E_0*DensityInc_i*SF*lambda_mean; + Jacobian_i[iVar][iVar] -= E_0*DensityInc_j*SF*lambda_mean; + } + + } + + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); + +} + +CAvgGradPBInc_Flow::CAvgGradPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + 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]; + } + + /*--- Primitive flow variables. ---*/ + + PrimVar_i = new su2double [nDim+4]; + PrimVar_j = new su2double [nDim+4]; + Mean_PrimVar = new su2double [nDim+4]; + + /*--- Incompressible flow, primitive variables nDim+2, (P, vx, vy, vz, rho) ---*/ + + Mean_GradPrimVar = new su2double*[nVar]; + + /*--- Incompressible flow, gradient primitive variables nDim+2, (P, vx, vy, vz, rho) ---*/ + + for (iVar = 0; iVar < nVar; iVar++) + Mean_GradPrimVar[iVar] = new su2double[nDim]; + +} + +CAvgGradPBInc_Flow::~CAvgGradPBInc_Flow(void) { + + delete [] PrimVar_i; + delete [] PrimVar_j; + delete [] Mean_PrimVar; + + for (iVar = 0; iVar < nVar; iVar++) + delete [] Mean_GradPrimVar[iVar]; + delete [] Mean_GradPrimVar; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + } + delete [] Jacobian_i; + delete [] Jacobian_j; + +} + +CNumerics::ResidualType<> CAvgGradPBInc_Flow::ComputeResidual(const CConfig *config) { + + /*--- Normalized normal vector ---*/ + + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = Normal[iDim]/Area; + + for (iVar = 0; iVar < nDim+4; iVar++) { + PrimVar_i[iVar] = V_i[iVar]; + PrimVar_j[iVar] = V_j[iVar]; + Mean_PrimVar[iVar] = 0.5*(PrimVar_i[iVar]+PrimVar_j[iVar]); + } + + /*--- Density and transport properties ---*/ + + Laminar_Viscosity_i = V_i[nDim+2]; Laminar_Viscosity_j = V_j[nDim+2]; + Eddy_Viscosity_i = V_i[nDim+3]; Eddy_Viscosity_j = V_j[nDim+3]; + + /*--- Mean transport properties ---*/ + + Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); + Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); + Mean_turb_ke = 0.0;//0.5*(turb_ke_i + turb_ke_j); + + /*--- Mean gradient approximation ---*/ + + for (iVar = 0; iVar < nVar; iVar++) + for (iDim = 0; iDim < nDim; iDim++) + Mean_GradPrimVar[iVar][iDim] = 0.5*(PrimVar_Grad_i[iVar+1][iDim] + PrimVar_Grad_j[iVar+1][iDim]); + + /*--- Get projected flux tensor ---*/ + + GetViscousPBIncProjFlux(Mean_PrimVar, Mean_GradPrimVar, Normal, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, Mean_turb_ke); + + + /*--- Implicit part ---*/ + + if (implicit) { + + dist_ij = 0.0; proj_vector_ij = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); + proj_vector_ij += (Coord_j[iDim]-Coord_i[iDim])*Normal[iDim]; + } + proj_vector_ij = proj_vector_ij/dist_ij; + dist_ij = sqrt(dist_ij); + + if (dist_ij == 0.0) { + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_j[iVar][jVar] = 0.0; + } + } + } + else { + GetViscousPBIncProjJacs(Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, dist_ij, UnitNormal, + Area, Jacobian_i, Jacobian_j); + } + + } + + return ResidualType<>(Proj_Flux_Tensor, Jacobian_i, Jacobian_j); +} + +CAvgGradCorrectedPBInc_Flow::CAvgGradCorrectedPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + 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]; + } + + /*--- Primitive flow variables. ---*/ + + Mean_PrimVar = new su2double [nDim+4]; + PrimVar_i = new su2double [nDim+4]; + PrimVar_j = new su2double [nDim+4]; + Proj_Mean_GradPrimVar_Edge = new su2double [nVar]; + Edge_Vector = new su2double [nDim]; + + Mean_GradPrimVar = new su2double* [nVar]; + for (iVar = 0; iVar < nVar; iVar++) + Mean_GradPrimVar[iVar] = new su2double [nDim]; + +} + +CAvgGradCorrectedPBInc_Flow::~CAvgGradCorrectedPBInc_Flow(void) { + + delete [] Mean_PrimVar; + delete [] PrimVar_i; + delete [] PrimVar_j; + delete [] Proj_Mean_GradPrimVar_Edge; + delete [] Edge_Vector; + + for (iVar = 0; iVar < nVar; iVar++) + delete [] Mean_GradPrimVar[iVar]; + delete [] Mean_GradPrimVar; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + } + delete [] Jacobian_i; + delete [] Jacobian_j; + +} + +CNumerics::ResidualType<> CAvgGradCorrectedPBInc_Flow::ComputeResidual(const CConfig *config) { + + unsigned short nPrimVarGrad = nDim+2, nPrimVar = nDim+4; + + /*--- Normalized normal vector ---*/ + + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = Normal[iDim]/Area; + + /*--- Conversion to Primitive Variables (P, u, v, w, rho, mu, muT) ---*/ + + for (iVar = 0; iVar < nPrimVar; iVar++) { + PrimVar_i[iVar] = V_i[iVar]; + PrimVar_j[iVar] = V_j[iVar]; + Mean_PrimVar[iVar] = 0.5*(PrimVar_i[iVar]+PrimVar_j[iVar]); + } + + /*--- Density and transport properties ---*/ + + DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; + Laminar_Viscosity_i = V_i[nDim+2]; Laminar_Viscosity_j = V_j[nDim+2]; + Eddy_Viscosity_i = V_i[nDim+3]; Eddy_Viscosity_j = V_j[nDim+3]; + + /*--- Mean transport properties ---*/ + + Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); + Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); + Mean_turb_ke = 0.0;//0.5*(turb_ke_i + turb_ke_j); + + /*--- Compute vector going from iPoint to jPoint ---*/ + + dist_ij_2 = 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]; + } + + /*--- Projection of the mean gradient in the direction of the edge ---*/ + + for (iVar = 0; iVar < nVar; iVar++) { + Proj_Mean_GradPrimVar_Edge[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPrimVar[iVar][iDim] = 0.5*(PrimVar_Grad_i[iVar+1][iDim] + PrimVar_Grad_j[iVar+1][iDim]); + Proj_Mean_GradPrimVar_Edge[iVar] += Mean_GradPrimVar[iVar][iDim]*Edge_Vector[iDim]; + } + if (dist_ij_2 != 0.0) { + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPrimVar[iVar][iDim] -= (Proj_Mean_GradPrimVar_Edge[iVar] - + (PrimVar_j[iVar+1]-PrimVar_i[iVar+1]))*Edge_Vector[iDim] / dist_ij_2; + } + } + } + + /*--- Get projected flux tensor ---*/ + + GetViscousPBIncProjFlux(Mean_PrimVar, Mean_GradPrimVar, Normal, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, Mean_turb_ke); + + /*--- Implicit part ---*/ + + if (implicit) { + + proj_vector_ij = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + proj_vector_ij += (Coord_j[iDim]-Coord_i[iDim])*Normal[iDim]; + } + proj_vector_ij = proj_vector_ij/dist_ij_2; + + if (dist_ij_2 == 0.0) { + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = 0.0; + Jacobian_j[iVar][jVar] = 0.0; + } + } + } + else { + GetViscousPBIncProjJacs(Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, sqrt(dist_ij_2), UnitNormal, + Area, Jacobian_i, Jacobian_j); + } + + } + + return ResidualType<>(Proj_Flux_Tensor, Jacobian_i, Jacobian_j); +} diff --git a/SU2_CFD/src/numerics/poisson.cpp b/SU2_CFD/src/numerics/poisson.cpp new file mode 100644 index 00000000000..f096a1b70fc --- /dev/null +++ b/SU2_CFD/src/numerics/poisson.cpp @@ -0,0 +1,421 @@ +/*! + * \file numerics_direct_poisson.cpp + * \brief This file contains all the convective term discretization. + * \author F. Palacios + * \version 6.1.0 "Falcon" + * + * The current SU2 release has been coordinated by the + * SU2 International Developers Society + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2018, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * 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/poisson.hpp" + + +CAvgGradCorrected_Poisson::CAvgGradCorrected_Poisson(unsigned short val_nDim, unsigned short val_nVar, + CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = false; + implicit = (config->GetKind_TimeIntScheme_Poisson() == EULER_IMPLICIT); + direct = false; + + + Edge_Vector = new su2double [nDim]; + Proj_Mean_GradPoissonVar_Edge = new su2double [nVar]; + Proj_Mean_GradPoissonVar_Kappa = new su2double [nVar]; + Proj_Mean_GradPoissonVar_Corrected = new su2double [nVar]; + Jacobian_i = new su2double* [nVar]; + Jacobian_j = new su2double* [nVar]; + Mean_GradPoissonVar = new su2double* [nVar]; + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nDim]; + Jacobian_j[iVar] = new su2double [nDim]; + Mean_GradPoissonVar[iVar] = new su2double [nDim]; + } + + Poisson_Coeff_i = 1.0; + Poisson_Coeff_j = 1.0; + + Mom_Coeff_i = new su2double [nDim]; + Mom_Coeff_j = new su2double [nDim]; + +} + + +CAvgGradCorrected_Poisson::~CAvgGradCorrected_Poisson(void) { + + delete [] Edge_Vector; + delete [] Proj_Mean_GradPoissonVar_Edge; + delete [] Proj_Mean_GradPoissonVar_Kappa; + delete [] Proj_Mean_GradPoissonVar_Corrected; + delete [] Mom_Coeff_i; + delete [] Mom_Coeff_j; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Mean_GradPoissonVar[iVar]; + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + } + delete [] Mean_GradPoissonVar; + delete [] Jacobian_i; + delete [] Jacobian_j; +} + +CNumerics::ResidualType<> CAvgGradCorrected_Poisson::ComputeResidual(const CConfig *config) { + + su2double Coeff_Mean; + su2double *MomCoeffxNormal = new su2double[nDim]; + su2double *MomCoeffxNormalCorrected = new su2double[nDim]; + su2double Mean_GradPoissonVar_Edge[MAXNDIM], GradPoisson[MAXNDIM]; + su2double Mean_GradPoissonVar_Face[MAXNDIM][MAXNDIM], Num, Den; + + Poisson_Coeff_Mean = 1.0;//0.5*(Poisson_Coeff_i + Poisson_Coeff_j); + + /*--- Multiply the normal with the coefficients of the momentum equation + * and use it instead of the normal during the projection operation ---*/ + for (iDim = 0; iDim < nDim; iDim++) { + Coeff_Mean = 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]) ; + MomCoeffxNormal[iDim] = Coeff_Mean*Normal[iDim]; + } + + /*--- Compute vector going from iPoint to jPoint ---*/ + /*--- Minimum correction approach. ---*/ + /*dist_ij_2 = 0; proj_vector_ij = 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]; + proj_vector_ij += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; + } + if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = proj_vector_ij/dist_ij_2;*/ + + /*--- Over relaxed approach. ---*/ + dist_ij_2 = 0; proj_vector_ij = 0; + Num = 0.0; Den = 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]; + Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; + Den += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; + } + if (Den == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = Num/Den; + + /*--- Orthogonal correction approach. ---*/ + /*dist_ij_2 = 0; proj_vector_ij = 0; + Num = 0.0; Den = 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]; + Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; + } + Den = dist_ij_2; + if (Den == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = sqrt(Num/Den);*/ + + /*--- Correct the face gradient for odd-even decoupling ---*/ + /*--- Steps are + * 1. Interpolate the gradient at the face -> du/ds|_in + * 2. Find the projection of the interpolated gradient on the edge vector -> (du/ds|_in) . e + * 3. Find the gradient as the difference between the neighboring nodes -> (u_j - u_i)/ds + * 4. Correct the gradient at the face using du/ds = du/ds|_in + ( (u_j - u_i)/ds - (du/ds|_in . e) ) . e ---*/ + + /*for (iVar = 0; iVar < nVar; iVar++) { + Mean_GradPoissonVar_Edge[iVar] = 0.0; + + GradPoisson[iVar] = (Poissonval_j - Poissonval_i)/sqrt(dist_ij_2); + + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); + + Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]/dist_ij_2; + } + } + + for (iVar = 0; iVar < nVar; iVar++) { + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPoissonVar_Face[iVar][iDim] = Mean_GradPoissonVar[iVar][iDim] + + (GradPoisson[iVar] - Mean_GradPoissonVar_Edge[iVar])*Edge_Vector[iDim]/dist_ij_2; + } + }*/ + + + /*--- Mean gradient approximation. Projection of the mean gradient + in the direction of the edge ---*/ + + for (iVar = 0; iVar < nVar; iVar++) { + Proj_Mean_GradPoissonVar_Edge[iVar] = 0.0; + Proj_Mean_GradPoissonVar_Kappa[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); + Proj_Mean_GradPoissonVar_Kappa[iVar] += Mean_GradPoissonVar[iVar][iDim]*MomCoeffxNormal[iDim]; + Proj_Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]; + } + Proj_Mean_GradPoissonVar_Corrected[iVar] = Proj_Mean_GradPoissonVar_Kappa[iVar]; + Proj_Mean_GradPoissonVar_Corrected[iVar] -= Proj_Mean_GradPoissonVar_Edge[iVar]*proj_vector_ij - + (Poissonval_j-Poissonval_i)*proj_vector_ij; + } + + /*--- Jacobians for implicit scheme ---*/ + + if (implicit) { + Jacobian_i[0][0] = -Poisson_Coeff_Mean*proj_vector_ij; + Jacobian_j[0][0] = Poisson_Coeff_Mean*proj_vector_ij; + } + + /*AD::SetPreaccOut(val_residual, nVar); + AD::EndPreacc();*/ + + delete [] MomCoeffxNormal; + delete [] MomCoeffxNormalCorrected; + + return ResidualType<>(Proj_Mean_GradPoissonVar_Corrected, Jacobian_i, Jacobian_j); +} + +CAvgGrad_Poisson::CAvgGrad_Poisson(unsigned short val_nDim, unsigned short val_nVar, + CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = false; + implicit = (config->GetKind_TimeIntScheme_Poisson() == EULER_IMPLICIT); // TODO: PBFlow: FIXED + direct = false; + + + Edge_Vector = new su2double [nDim]; + Proj_Mean_GradPoissonVar_Normal = new su2double [nVar]; + Proj_Mean_GradPoissonVar_Corrected = new su2double [nVar]; + Mom_Coeff_i = new su2double [nDim]; + Mom_Coeff_j = new su2double [nDim]; + + Jacobian_i = new su2double* [nVar]; + Jacobian_j = new su2double* [nVar]; + Mean_GradPoissonVar = new su2double* [nVar]; + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nDim]; + Jacobian_j[iVar] = new su2double [nDim]; + Mean_GradPoissonVar[iVar] = new su2double [nDim]; + } + + Poisson_Coeff_i = 1.0; + Poisson_Coeff_j = 1.0; +} + + +CAvgGrad_Poisson::~CAvgGrad_Poisson(void) { + + unsigned int iDim; + delete [] Edge_Vector; + delete [] Proj_Mean_GradPoissonVar_Normal; + delete [] Proj_Mean_GradPoissonVar_Corrected; + delete [] Mom_Coeff_i; + delete [] Mom_Coeff_j; + + for (iVar = 0; iVar < nVar; iVar++) { + delete [] Mean_GradPoissonVar[iVar]; + delete [] Jacobian_i[iVar]; + delete [] Jacobian_j[iVar]; + } + delete [] Mean_GradPoissonVar; + delete [] Jacobian_i; + delete [] Jacobian_j; + +} + +CNumerics::ResidualType<> CAvgGrad_Poisson::ComputeResidual(const CConfig *config) { + + + su2double Coeff_Mean, Num, Den; + su2double *MomCoeffxNormal = new su2double[nDim]; + su2double Mean_GradPoissonVar_Edge[MAXNDIM], GradPoisson[MAXNDIM]; + su2double Mean_GradPoissonVar_Face[MAXNDIM][MAXNDIM]; + + + + Poisson_Coeff_Mean = 1.0;//0.5*(Poisson_Coeff_i + Poisson_Coeff_j); + + /*--- Multiply the normal with the coefficients of the momentum equation + * and use it instead of the normal during the projection operation ---*/ + for (iDim = 0; iDim < nDim; iDim++) { + Coeff_Mean = 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]) ; + MomCoeffxNormal[iDim] = Coeff_Mean*Normal[iDim]; + } + + /*--- Compute vector going from iPoint to jPoint ---*/ + /*--- Minimum correction approach. ---*/ + /*dist_ij_2 = 0; proj_vector_ij = 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]; + proj_vector_ij += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; + } + if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = proj_vector_ij/dist_ij_2;*/ + + /*--- Over relaxed approach. ---*/ + dist_ij_2 = 0; proj_vector_ij = 0; + Num = 0.0; Den = 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]; + Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; + Den += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; + } + if (Den == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = Num/Den; + + /*--- Orthogonal correction approach. ---*/ + /*dist_ij_2 = 0; proj_vector_ij = 0; + Num = 0.0; Den = 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]; + Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; + } + Den = dist_ij_2; + if (Den == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = sqrt(Num/Den);*/ + + /*--- Correct the face gradient for odd-even decoupling ---*/ + /*--- Steps are + * 1. Interpolate the gradient at the face -> du/ds|_in + * 2. Find the projection of the interpolated gradient on the edge vector -> (du/ds|_in) . e + * 3. Find the gradient as the difference between the neighboring nodes -> (u_j - u_i)/ds + * 4. Correct the gradient at the face using du/ds = du/ds|_in + ( (u_j - u_i)/ds - (du/ds|_in . e) ) . e ---*/ + + for (iVar = 0; iVar < nVar; iVar++) { + Mean_GradPoissonVar_Edge[iVar] = 0.0; + + GradPoisson[iVar] = (Poissonval_j - Poissonval_i)/sqrt(dist_ij_2); + + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); + + Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]/sqrt(dist_ij_2); + } + } + + for (iVar = 0; iVar < nVar; iVar++) { + for (iDim = 0; iDim < nDim; iDim++) { + Mean_GradPoissonVar_Face[iVar][iDim] = Mean_GradPoissonVar[iVar][iDim] + + (GradPoisson[iVar] - Mean_GradPoissonVar_Edge[iVar])*Edge_Vector[iDim]/sqrt(dist_ij_2); + } + } + + /*--- Mean gradient approximation. Projection of the mean gradient in the direction of the edge ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + Proj_Mean_GradPoissonVar_Normal[iVar] = 0.0; + Proj_Mean_GradPoissonVar_Corrected[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Proj_Mean_GradPoissonVar_Normal[iVar] += Mean_GradPoissonVar_Face[iVar][iDim]*MomCoeffxNormal[iDim]; + } + } + + /*--- Jacobians for implicit scheme ---*/ + if (implicit) { + Jacobian_i[0][0] = -Poisson_Coeff_Mean*proj_vector_ij; + Jacobian_j[0][0] = Poisson_Coeff_Mean*proj_vector_ij; + } + + delete [] MomCoeffxNormal; + + return ResidualType<>(Proj_Mean_GradPoissonVar_Normal, Jacobian_i, Jacobian_j); +} + + + +CPressure_Poisson::CPressure_Poisson(unsigned short val_nDim, unsigned short val_nVar, + CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = true; + + Edge_Vector = new su2double [nDim]; + + Poisson_Coeff_i = 1.0; + Poisson_Coeff_j = 1.0; + + Mom_Coeff_i = new su2double [nDim]; + Mom_Coeff_j = new su2double [nDim]; + +} + + +CPressure_Poisson::~CPressure_Poisson(void) { + + unsigned int iDim; + delete [] Edge_Vector; + + delete Mom_Coeff_i; + delete Mom_Coeff_j; +} + +void CPressure_Poisson::ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) { + + + su2double Area; + su2double UnitNormal[MAXNDIM]; + + /*--- Compute vector going from iPoint to jPoint ---*/ + + dist_ij_2 = 0; proj_vector_ij = 0; + Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + for (iDim = 0; iDim < nDim; iDim++) { + Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; + dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; + proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; + UnitNormal[iDim] = Normal[iDim]/Area; + } + if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; + else proj_vector_ij = proj_vector_ij/dist_ij_2; + + + Poisson_Coeff_Mean = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Poisson_Coeff_Mean += 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]); + + + val_residual[0] = 0.5*Poisson_Coeff_Mean*(Poissonval_i + Poissonval_j); + + if (implicit){ + Jacobian_i[0][0] = Poisson_Coeff_Mean; + Jacobian_j[0][0] = -Poisson_Coeff_Mean; + } +} + + +CSource_PoissonFVM::CSource_PoissonFVM(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { } + +CSource_PoissonFVM::~CSource_PoissonFVM(void) { } + +void CSource_PoissonFVM::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { + + if (config->GetKind_Incomp_System()==ENUM_INCOMP_SYSTEM::PRESSURE_BASED) + val_residual[0] = Source_Term; + else + val_residual[0] = Source_Term*Volume; +} From 8d0bcf91c4c159ed196c8a9b04b37ffd0e70f11c Mon Sep 17 00:00:00 2001 From: Thijs Aalbers Date: Tue, 12 May 2026 13:24:43 +0200 Subject: [PATCH 3/4] Remove outdated numerics --- SU2_CFD/include/numerics/CNumerics.hpp | 103 ----- SU2_CFD/include/numerics/pbflow.hpp | 232 ---------- SU2_CFD/include/numerics/poisson.hpp | 219 ---------- SU2_CFD/src/drivers/CDriver.cpp | 36 +- SU2_CFD/src/meson.build | 2 - SU2_CFD/src/numerics/CNumerics.cpp | 195 +-------- SU2_CFD/src/numerics/pbflow.cpp | 573 ------------------------- SU2_CFD/src/numerics/poisson.cpp | 421 ------------------ 8 files changed, 18 insertions(+), 1763 deletions(-) delete mode 100644 SU2_CFD/include/numerics/pbflow.hpp delete mode 100644 SU2_CFD/include/numerics/poisson.hpp delete mode 100644 SU2_CFD/src/numerics/pbflow.cpp delete mode 100644 SU2_CFD/src/numerics/poisson.cpp diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 9547164d92a..535e40cba96 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -57,7 +57,6 @@ class CNumerics { su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl's number. */ su2double MassFlux; /*!< \brief Mass flux across edge. */ su2double - **Flux_Tensor, /*!< \brief Flux tensor (used for viscous and inviscid purposes). */ *Proj_Flux_Tensor; /*!< \brief Flux tensor projected in a direction. */ su2double **tau; /*!< \brief Viscous stress tensor. */ const su2double delta [3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; /*!< \brief Identity matrix. */ @@ -157,12 +156,6 @@ class CNumerics { TurbPsi_Grad_j, /*!< \brief Gradient of adjoint turbulent variables at point j. */ AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ - su2double - Poisson_Coeff_i, /*!< \brief Poisson coefficient at point i. */ - Poisson_Coeff_j; /*!< \brief Poisson coefficient at point j. */ - su2double - Poissonval_i, /*!< \brief Poisson value at point i. */ - Poissonval_j; /*!< \brief Poisson value at point j. */ su2double LocalGridLength_i; /*!< \brief Local grid length at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ @@ -187,7 +180,6 @@ class CNumerics { su2double StrainMag_i, StrainMag_j; /*!< \brief Strain rate magnitude. */ su2double Dissipation_i, Dissipation_j; /*!< \brief Dissipation. */ su2double Dissipation_ij; - su2double Source_Term; /*!< \brief Source term for poisson equation. */ su2double roughness_i = 0.0, /*!< \brief Roughness of the wall nearest to point i. */ roughness_j = 0.0; /*!< \brief Roughness of the wall nearest to point j. */ @@ -207,7 +199,6 @@ class CNumerics { su2double uq_delta_b; /*!< \brief Magnitude of perturbation */ su2double uq_urlx; /*!< \brief Under-relaxation factor for numerical stability */ bool uq_permute; /*!< \brief Flag for eigenvector permutation */ - su2double Face_Flux; /*!< \brief Face flux */ bool nemo; /*!< \brief Flag for NEMO problems */ @@ -1113,40 +1104,6 @@ class CNumerics { */ inline su2double GetDissipation() const { return Dissipation_ij; } - /*! - * \brief Set the Poisson coefficient - * \param[in] val_Poisson_i - Value of the Poisson coefficient at point i. - * \param[in] val_Poisson_j - Value of the Poisson coefficient at point j. - */ - inline void SetPoisson_Coeff(su2double val_Poisson_Coeff_i,su2double val_Poisson_Coeff_j){ - Poisson_Coeff_i = val_Poisson_Coeff_i; - Poisson_Coeff_j = val_Poisson_Coeff_j; - } - - /*! - * \brief Set the Poisson value - * \param[in] val_Poisson_i - Value of the Poisson variable at point i. - * \param[in] val_Poisson_j - Value of the Poisson variable at point j. - */ - inline void SetPoissonval(su2double val_Poisson_i,su2double val_Poisson_j){ - Poissonval_i = val_Poisson_i; - Poissonval_j = val_Poisson_j; - } - - /*! - * \brief Set the Poisson value - * \param[in] val_Poisson_i - Value of the Poisson variable at point i. - * \param[in] val_Poisson_j - Value of the Poisson variable at point j. - */ - inline void SetSourcePoisson(su2double val_Source_Term) { Source_Term = val_Source_Term; } - - /*! - * \brief Set the value of the momentum equation coefficients for Poisson eq. - * \param[in] val_Mom_Coeff_i - Value of the cross coefficient at point i. - * \param[in] val_Mom_Coeff_j - Value of the cross coefficient at point j. - */ - inline virtual void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { } - /*! * \brief Compute the projected inviscid flux vector. * \param[in] val_density - Pointer to the density. @@ -1220,21 +1177,6 @@ class CNumerics { su2double val_scale, su2double **val_Proj_Jac_Tensor) const; - /*! - * \brief Compute the projected inviscid flux vector for (pressure based) incompresible simulations - * \param[in] val_density - Pointer to the density. - * \param[in] val_velocity - Pointer to the velocity. - * \param[in] val_pressure - Pointer to the pressure. - * \param[in] val_betainc2 - Value of the artificial compresibility factor. - * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. - * \param[out] val_Proj_Flux - Pointer to the projected flux. - */ - void GetInviscidPBProjFlux(const su2double *val_density, - const su2double *val_velocity, - const su2double *val_pressure, - const su2double *val_normal, - su2double *val_Proj_Flux); - /*! * \brief Compute the low speed preconditioning matrix. * \param[in] val_density - Value of the density. @@ -1288,19 +1230,6 @@ class CNumerics { const su2double *val_Mean_SecVar, su2double **val_Jac_PC) const; - /*! - * \brief Compute the projection of the inviscid Jacobian matrices (Pressure based method). - * \param[in] val_density - Value of the density. - * \param[in] val_velocity - Pointer to the velocity. - * \param[in] val_cp - Value of the specific heat at constant pressure. - * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. - * \param[in] val_scale - Scale of the projection. - * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. - */ - void GetInviscidPBProjJac(const su2double val_density, const su2double *val_velocity, - const su2double *val_normal, const su2double val_scale, - su2double **val_Proj_Jac_tensor); - /*! * \overload * \brief Computation of the matrix P for a generic fluid model @@ -1925,38 +1854,6 @@ class CNumerics { * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport */ inline bool GetBoundedScalar() const { return bounded_scalar;} - - /*! - * \brief Compute the projection of the viscous Jacobian matrices. - * \param[in] val_laminar_viscosity - Value of the laminar viscosity. - * \param[in] val_eddy_viscosity - Value of the eddy viscosity. - * \param[in] val_dist_ij - Distance between the points. - * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. - * \param[in] val_dS - Area of the face between two nodes. - * \param[out] val_Proj_Jac_Tensor_i - Pointer to the projected viscous Jacobian at point i. - * \param[out] val_Proj_Jac_Tensor_j - Pointer to the projected viscous Jacobian at point j. - */ - void GetViscousPBIncProjJacs(const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity, - const su2double val_dist_ij, - const su2double *val_normal, - const su2double val_dS, - su2double **val_Proj_Jac_Tensor_i, - su2double **val_Proj_Jac_Tensor_j); - /*! - * \brief Compute the projection of the viscous fluxes into a direction (pressure based method). - * \param[in] val_primvar - Primitive variables. - * \param[in] val_gradprimvar - Gradient of the primitive variables. - * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. - * \param[in] val_laminar_viscosity - Laminar viscosity. - * \param[in] val_eddy_viscosity - Eddy viscosity. - */ - void GetViscousPBIncProjFlux(const su2double *val_primvar, - su2double **val_gradprimvar, - const su2double *val_normal, - const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity, - const su2double val_turb_ke); }; /*! diff --git a/SU2_CFD/include/numerics/pbflow.hpp b/SU2_CFD/include/numerics/pbflow.hpp deleted file mode 100644 index 637db426d1f..00000000000 --- a/SU2_CFD/include/numerics/pbflow.hpp +++ /dev/null @@ -1,232 +0,0 @@ -/*! - * \file pbflow.hpp - * \brief Delarations of numerics classes for heat transfer problems. - * \author F. Palacios, T. Economon - * \version 8.0.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2020, 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 CUpwPB_Flow - * \brief Class for solving an approximate Riemann solver of Roe for the pressure based incompressible flow equations. - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CUpwPB_Flow : public CNumerics { -private: - bool implicit, dynamic_grid; - bool gravity; - su2double Froude, Upw_i, Upw_j; - su2double *Diff_U; - su2double *Velocity_i, *Velocity_j, *MeanVelocity, *Velocity_upw; - su2double *ProjFlux_i, *ProjFlux_j; - su2double Proj_ModJac_Tensor_ij, Pressure_i, - Pressure_j, MeanDensity, MeanSoundSpeed, MeanPressure, MeanBetaInc2, - ProjVelocity, FaceVel, Face_Flux; - unsigned short iDim, iVar, jVar, kVar; - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - su2double **Jacobian_upw = nullptr; - -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CUpwPB_Flow(void); - - /*! - * \brief Compute the Roe's flux between two nodes i and j. - * \param[out] val_residual - Pointer to the total 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) override; - - /*! - * \brief Set the value of face velocity. This is used as a proxy for massflux at a face. - * \param[in] val_FaceVel. - */ - inline void SetFaceVel(su2double val_FaceVel){ FaceVel = val_FaceVel; } -}; - - -/*! - * \class CCentLaxArtComp_Flow - * \brief Class for computing the Lax-Friedrich centered scheme (artificial compressibility). - * \ingroup ConvDiscr - * \author F. Palacios - */ -class CCentPB_Flow : public CNumerics { -private: - bool implicit, dynamic_grid; - bool gravity; - su2double Froude, Upw_i, Upw_j; - su2double *Diff_U; - su2double *Velocity_i, *Velocity_j, *MeanVelocity, *Velocity_upw; - su2double *ProjFlux_i, *ProjFlux_j; - su2double *Lambda, *Epsilon; - su2double **P_Tensor, **invP_Tensor,**val_Jacobian_upw; - su2double Proj_ModJac_Tensor_ij, Pressure_i, - Pressure_j, MeanDensity, MeanSoundSpeed, MeanPressure, MeanBetaInc2, - ProjVelocity, FaceVel; - unsigned short iDim, iVar, jVar, kVar; - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - su2double **Jacobian_upw = nullptr; - -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. - */ - CCentPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CCentPB_Flow(void); - - /*! - * \brief Compute the flow residual using a Lax method. - * \param[out] val_resconv - Pointer to the convective residual. - * \param[out] val_resvisc - Pointer to the artificial viscosity 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) override; -}; - -/*! - * \class CAvgGradInc_Flow - * \brief Class for computing viscous term using an average of gradients. - * \ingroup ViscDiscr - * \author A. Bueno, F. Palacios, T. Economon - */ -class CAvgGradPBInc_Flow : public CNumerics { -private: - unsigned short iDim, iVar, jVar; /*!< \brief Iterators in dimension an variable. */ - su2double *Mean_PrimVar, /*!< \brief Mean primitive variables. */ - *PrimVar_i, *PrimVar_j; /*!< \brief Primitives variables at point i and j. */ - su2double **Mean_GradPrimVar, /*!< \brief Mean value of the gradient. */ - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, /*!< \brief Mean value of the viscosity. */ - Mean_turb_ke, /*!< \brief Mean value of the turbulent kinetic energy. */ - dist_ij, /*!< \brief Length of the edge and face. */ - proj_vector_ij; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - bool implicit; /*!< \brief Implicit calculus. */ - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - -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. - */ - CAvgGradPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGradPBInc_Flow(void); - /*! - * \brief Compute the viscous flow residual using an average of gradients. - * \param[out] val_residual - Pointer to the total 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) override; -}; - - -/*! - * \class CAvgGradCorrectedPBInc_Flow - * \brief Class for computing viscous term using an average of gradients with correction (incompressible). - * \ingroup ViscDiscr - */ -class CAvgGradCorrectedPBInc_Flow : public CNumerics { -private: - unsigned short iDim, iVar, jVar; /*!< \brief Iterators in dimension an variable. */ - su2double *Mean_PrimVar; /*!< \brief Mean primitive variables. */ - su2double *PrimVar_i, *PrimVar_j, /*!< \brief Primitives variables at point i and 1. */ - *Edge_Vector, /*!< \brief Vector form point i to point j. */ - **Mean_GradPrimVar, *Proj_Mean_GradPrimVar_Edge, /*!< \brief Mean value of the gradient. */ - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, /*!< \brief Mean value of the viscosity. */ - Mean_turb_ke, /*!< \brief Mean value of the turbulent kinetic energy. */ - dist_ij_2, /*!< \brief Length of the edge and face. */ - proj_vector_ij; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - bool implicit; /*!< \brief Implicit calculus. */ - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - -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. - */ - CAvgGradCorrectedPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGradCorrectedPBInc_Flow(void); - - /*! - * \brief Compute the viscous flow residual using an average of gradients with correction. - * \param[out] val_residual - Pointer to the total 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) override; -}; diff --git a/SU2_CFD/include/numerics/poisson.hpp b/SU2_CFD/include/numerics/poisson.hpp deleted file mode 100644 index 21fe28d60b0..00000000000 --- a/SU2_CFD/include/numerics/poisson.hpp +++ /dev/null @@ -1,219 +0,0 @@ -/*! - * \file heat.hpp - * \brief Delarations of numerics classes for heat transfer problems. - * \author F. Palacios, T. Economon - * \version 8.0.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2020, 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 CAvgGrad_Poisson - * \brief Class for computing viscous term using average of gradients without correction (Poisson equation). - * \ingroup ViscDiscr - */ -class CAvgGrad_Poisson : public CNumerics { -private: - - su2double *Edge_Vector; - bool implicit,direct; - su2double **Mean_GradPoissonVar; - su2double *Mom_Coeff_i,*Mom_Coeff_j; - su2double *Proj_Mean_GradPoissonVar_Normal, *Proj_Mean_GradPoissonVar_Corrected; - su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean ; - unsigned short iVar, iDim; - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGrad_Poisson(void); - - /*! - * \brief Compute the viscous heat residual using an average of gradients with correction. - * \param[out] val_residual - Pointer to the total residual. - * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[out] 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) override; - - - inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; - Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; - } - } -}; - -/*! - * \class CAvgGradCorrected_Poisson - * \brief Class for computing viscous term using average of gradients with correction (Poisson equation). - * \ingroup ViscDiscr - */ - -class CAvgGradCorrected_Poisson : public CNumerics { -private: - - su2double *Edge_Vector; - bool implicit,direct; - su2double **Mean_GradPoissonVar; - su2double *Mom_Coeff_i,*Mom_Coeff_j; - su2double *Proj_Mean_GradPoissonVar_Kappa, *Proj_Mean_GradPoissonVar_Edge, *Proj_Mean_GradPoissonVar_Corrected; - su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean; - unsigned short iVar, iDim; - - su2double *Flux = nullptr; - su2double **Jacobian_i = nullptr; - su2double **Jacobian_j = nullptr; - -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CAvgGradCorrected_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CAvgGradCorrected_Poisson(void); - - /*! - * \brief Compute the viscous residual using an average of gradients with correction. - * \param[out] val_residual - Pointer to the total residual. - * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[out] 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) override; - - - inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; - Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; - } - } -}; - -/*! - * \class CAvgGrad_Poisson - * \brief Class for computing viscous term using average of gradients without correction (Poisson equation). - * \ingroup ViscDiscr - */ -class CPressure_Poisson : public CNumerics { -private: - - su2double *Edge_Vector; - bool implicit,direct; - su2double **Mean_GradPoissonVar; - su2double *Proj_Mean_GradPoissonVar_Kappa, *Proj_Mean_GradPoissonVar_Edge, *Proj_Mean_GradPoissonVar_Corrected; - su2double dist_ij_2, proj_vector_ij, Poisson_Coeff_Mean; - su2double *Mom_Coeff_i,*Mom_Coeff_j; - unsigned short iVar, iDim; - -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CPressure_Poisson(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CPressure_Poisson(void); - - /*! - * \brief Compute the viscous heat residual using an average of gradients with correction. - * \param[out] val_residual - Pointer to the total residual. - * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[out] Jacobian_j - Jacobian of the numerical method at node j (implicit computation). - * \param[in] config - Definition of the particular problem. - */ - void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config); - - - inline void SetInvMomCoeff(su2double *val_Mom_Coeff_i, su2double *val_Mom_Coeff_j) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - Mom_Coeff_i[iDim] = val_Mom_Coeff_i[iDim]; - Mom_Coeff_j[iDim] = val_Mom_Coeff_j[iDim]; - } - } -}; - -/*! - * \class CSource_Poisson - * \brief Class for source term of the Poisson equation. - * \ingroup SourceDiscr - */ -class CSource_PoissonFVM : public CNumerics { -public: - - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Name of the input config file - * - */ - CSource_PoissonFVM(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); - - - /*! - * \brief Residual for source term integration. - * \param[out] val_residual - Pointer to the total residual. - * \param[out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). - * \param[in] config - Definition of the particular problem. - */ - void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, CConfig *config); - - /*! - * \brief Destructor of the class. - */ - ~CSource_PoissonFVM(void); -}; - diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index c95474c9ef7..b2f4ec1360a 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -56,8 +56,6 @@ #include "../../include/numerics/template.hpp" #include "../../include/numerics/radiation.hpp" #include "../../include/numerics/heat.hpp" -#include "../../include/numerics/pbflow.hpp" -#include "../../include/numerics/poisson.hpp" #include "../../include/numerics/flow/convection/roe.hpp" #include "../../include/numerics/flow/convection/fds.hpp" #include "../../include/numerics/flow/convection/fvs.hpp" @@ -1680,16 +1678,16 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver } else { /*--- Incompressible flow, use pressure-based method ---*/ switch (config->GetKind_Centered_Flow()) { - case CENTERED::CDS : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentPB_Flow(nDim, nVar_Flow, config); break; + case CENTERED::CDS : /* numerics[MESH_0][FLOW_SOL][conv_term] = new CCentPB_Flow(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 CCentPB_Flow(nDim, nVar_Flow, config); + // for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + // numerics[iMGlevel][FLOW_SOL][conv_term] = new CCentPB_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 CUpwPB_Flow(nDim, nVar_Flow, config); + // for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + // numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwPB_Flow(nDim, nVar_Flow, config); } } @@ -1816,8 +1814,8 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver switch (config->GetKind_Upwind_Flow()) { case UPWIND::UDS: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwPB_Flow(nDim, nVar_Flow, config); - numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwPB_Flow(nDim, nVar_Flow, config); + // numerics[iMGlevel][FLOW_SOL][conv_term] = new CUpwPB_Flow(nDim, nVar_Flow, config); + // numerics[iMGlevel][FLOW_SOL][conv_bound_term] = new CUpwPB_Flow(nDim, nVar_Flow, config); } break; default: @@ -1871,13 +1869,13 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradInc_Flow(nDim, nVar_Flow, false, config); } else { /*--- Incompressible flow, use pressure based method ---*/ - numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradCorrectedPBInc_Flow(nDim, nVar_Flow, config); - for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); + // numerics[MESH_0][FLOW_SOL][visc_term] = new CAvgGradCorrectedPBInc_Flow(nDim, nVar_Flow, config); + // for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + // numerics[iMGlevel][FLOW_SOL][visc_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); /*--- Definition of the boundary condition method ---*/ - for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); + // for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + // numerics[iMGlevel][FLOW_SOL][visc_bound_term] = new CAvgGradPBInc_Flow(nDim, nVar_Flow, config); } } @@ -2150,18 +2148,18 @@ 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 CAvgGradCorrected_Poisson(nDim, nVar_Poisson, config); + // numerics[MESH_0][POISSON_SOL][visc_term] = new CAvgGradCorrected_Poisson(nDim, nVar_Poisson, config); for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) - numerics[iMGlevel][POISSON_SOL][visc_term] = new CAvgGrad_Poisson(nDim, nVar_Poisson, config); + // numerics[iMGlevel][POISSON_SOL][visc_term] = new CAvgGrad_Poisson(nDim, nVar_Poisson, config); /*--- 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_Poisson(nDim, nVar_Poisson, config); + // numerics[iMGlevel][POISSON_SOL][visc_bound_term] = new CAvgGrad_Poisson(nDim, nVar_Poisson, config); /*--- Definition of the source term integration scheme for each equation and mesh level ---*/ - numerics[iMGlevel][POISSON_SOL][source_first_term] = new CSource_PoissonFVM(nDim, nVar_Poisson, config); - numerics[iMGlevel][POISSON_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Poisson, config); + // numerics[iMGlevel][POISSON_SOL][source_first_term] = new CSource_PoissonFVM(nDim, nVar_Poisson, config); + // numerics[iMGlevel][POISSON_SOL][source_second_term] = new CSourceNothing(nDim, nVar_Poisson, config); } } diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 6b754aec007..3b40822a34f 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -122,8 +122,6 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', su2_cfd_src += files(['numerics/CNumerics.cpp', 'numerics/template.cpp', 'numerics/radiation.cpp', - 'numerics/poisson.cpp', - 'numerics/pbflow.cpp', 'numerics/flow/convection/roe.cpp', 'numerics/flow/convection/fds.cpp', 'numerics/flow/convection/fvs.cpp', diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp index 6ff3efcc9aa..bd8728e0e83 100644 --- a/SU2_CFD/src/numerics/CNumerics.cpp +++ b/SU2_CFD/src/numerics/CNumerics.cpp @@ -44,7 +44,7 @@ CNumerics::CNumerics() { CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) { - unsigned short iDim, iVar; + unsigned short iDim; nDim = val_nDim; nVar = val_nVar; @@ -55,10 +55,6 @@ CNumerics::CNumerics(unsigned short val_nDim, unsigned short val_nVar, Gas_Constant = config->GetGas_ConstantND(); energy_multicomponent = config->GetKind_FluidModel() == FLUID_MIXTURE && config->GetEnergy_Equation(); - Flux_Tensor = new su2double* [nVar]; - for (iVar = 0; iVar < (nVar); iVar++) - Flux_Tensor[iVar] = new su2double [nDim] (); - tau = new su2double* [nDim]; for (iDim = 0; iDim < nDim; iDim++) tau[iDim] = new su2double [nDim] (); @@ -84,12 +80,6 @@ CNumerics::~CNumerics() { // visc delete [] Proj_Flux_Tensor; - if (Flux_Tensor) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) - delete [] Flux_Tensor[iVar]; - delete [] Flux_Tensor; - } - if (tau) { for (unsigned short iDim = 0; iDim < nDim; iDim++) delete [] tau[iDim]; @@ -1585,186 +1575,3 @@ su2double CNumerics::GetRoe_Dissipation(const su2double Dissipation_i, } return Dissipation_ij; } - -void CNumerics::GetInviscidPBProjFlux(const su2double *val_density, - const su2double *val_velocity, - const su2double *val_pressure, - const su2double *val_normal, - su2double *val_Proj_Flux) { - su2double rhou, rhov, rhow; - - if (nDim == 2) { - rhou = (*val_density)*val_velocity[0]; - rhov = (*val_density)*val_velocity[1]; - - val_Proj_Flux[0] = (rhou*val_velocity[0])*val_normal[0] + rhou*val_velocity[1]*val_normal[1]; - val_Proj_Flux[1] = rhov*val_velocity[0]*val_normal[0] + (rhov*val_velocity[1])*val_normal[1]; - } - else { - rhou = (*val_density)*val_velocity[0]; - rhov = (*val_density)*val_velocity[1]; - rhow = (*val_density)*val_velocity[2]; - - val_Proj_Flux[0] = (rhou*val_velocity[0])*val_normal[0] + rhou*val_velocity[1]*val_normal[1] + rhou*val_velocity[2]*val_normal[2]; - val_Proj_Flux[1] = rhov*val_velocity[0]*val_normal[0] + (rhov*val_velocity[1])*val_normal[1] + rhov*val_velocity[2]*val_normal[2]; - val_Proj_Flux[2] = rhow*val_velocity[0]*val_normal[0] + rhow*val_velocity[1]*val_normal[1] + (rhow*val_velocity[2])*val_normal[2]; - } -} - -void CNumerics::GetInviscidPBProjJac(const su2double val_density, const su2double *val_velocity, const su2double *val_normal, - const su2double val_scale, su2double **val_Proj_Jac_Tensor) { - - unsigned short iDim; - su2double proj_vel; - - proj_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - proj_vel += val_velocity[iDim]*val_normal[iDim]; - - if (nDim == 2) { - val_Proj_Jac_Tensor[0][0] = val_scale*val_density*(val_velocity[0]*val_normal[0] + proj_vel); - val_Proj_Jac_Tensor[0][1] = val_scale*val_density*val_velocity[0]*val_normal[1]; - - val_Proj_Jac_Tensor[1][0] = val_scale*val_density*val_velocity[1]*val_normal[0]; - val_Proj_Jac_Tensor[1][1] = val_scale*val_density*(val_velocity[1]*val_normal[1] + proj_vel); - } - else { - val_Proj_Jac_Tensor[0][0] = val_scale*val_density*(proj_vel+val_velocity[0]*val_normal[0]); - val_Proj_Jac_Tensor[0][1] = val_scale*val_density*(val_velocity[0]*val_normal[1]); - val_Proj_Jac_Tensor[0][2] = val_scale*val_density*(val_velocity[0]*val_normal[2]); - - val_Proj_Jac_Tensor[1][0] = val_scale*val_density*(val_velocity[1]*val_normal[0]); - val_Proj_Jac_Tensor[1][1] = val_scale*val_density*(proj_vel+val_velocity[1]*val_normal[1]); - val_Proj_Jac_Tensor[1][2] = val_scale*val_density*(val_velocity[1]*val_normal[2]); - - val_Proj_Jac_Tensor[2][0] = val_scale*val_density*(val_velocity[2]*val_normal[0]); - val_Proj_Jac_Tensor[2][1] = val_scale*val_density*(val_velocity[2]*val_normal[1]); - val_Proj_Jac_Tensor[2][2] = val_scale*val_density*(proj_vel+val_velocity[2]*val_normal[2]); - - } -} - -void CNumerics::GetViscousPBIncProjFlux(const su2double *val_primvar, - su2double **val_gradprimvar, - const su2double *val_normal, - const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity, - const su2double val_turb_ke) { - - unsigned short iVar, iDim, jDim; - su2double total_viscosity, div_vel, Density; - - Density = val_primvar[nDim+1]; - - total_viscosity = (val_laminar_viscosity + val_eddy_viscosity); - - /*--- The full stress tensor is needed for variable density, as nabla.u != 0 ---*/ - /*--- Note: Gradients are computed as [iDim][jDim] and not [iDim+1][jDim] because - *--- the mean gradient passed here contains only velocities and no pressure. The - *--- mean gradient is computed as PrimVar_Grad[iVar+1][iDim], so no need to repeat.---*/ - - div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - div_vel += val_gradprimvar[iDim][iDim]; - - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = (total_viscosity*(val_gradprimvar[jDim][iDim] + - val_gradprimvar[iDim][jDim] ) - -TWO3*total_viscosity*div_vel*delta[iDim][jDim] - -TWO3*Density*val_turb_ke*delta[iDim][jDim]); - - /*--- Gradient of primitive variables -> [Pressure vel_x vel_y vel_z Density] ---*/ - - if (nDim == 2) { - Flux_Tensor[0][0] = tau[0][0]; - Flux_Tensor[1][0] = tau[0][1]; - - Flux_Tensor[0][1] = tau[1][0]; - Flux_Tensor[1][1] = tau[1][1]; - - } else { - - Flux_Tensor[0][0] = tau[0][0]; - Flux_Tensor[1][0] = tau[0][1]; - Flux_Tensor[2][0] = tau[0][2]; - - Flux_Tensor[0][1] = tau[1][0]; - Flux_Tensor[1][1] = tau[1][1]; - Flux_Tensor[2][1] = tau[1][2]; - - Flux_Tensor[0][2] = tau[2][0]; - Flux_Tensor[1][2] = tau[2][1]; - Flux_Tensor[2][2] = tau[2][2]; - - } - - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Flux_Tensor[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Proj_Flux_Tensor[iVar] += Flux_Tensor[iVar][iDim] * val_normal[iDim]; - } - -} - -void CNumerics::GetViscousPBIncProjJacs(const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity, - const su2double val_dist_ij, - const su2double *val_normal, - const su2double val_dS, - su2double **val_Proj_Jac_Tensor_i, su2double **val_Proj_Jac_Tensor_j) { - unsigned short iDim, iVar, jVar; - - su2double theta = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - theta += val_normal[iDim]*val_normal[iDim]; - - su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity; - su2double factor = total_viscosity/(val_dist_ij)*val_dS; - - if (nDim == 3) { - su2double thetax = theta + val_normal[0]*val_normal[0]/3.0; - su2double thetay = theta + val_normal[1]*val_normal[1]/3.0; - su2double thetaz = theta + val_normal[2]*val_normal[2]/3.0; - - su2double etax = val_normal[1]*val_normal[2]/3.0; - su2double etay = val_normal[0]*val_normal[2]/3.0; - su2double etaz = val_normal[0]*val_normal[1]/3.0; - - val_Proj_Jac_Tensor_i[0][0] = -factor*thetax; - val_Proj_Jac_Tensor_i[0][1] = -factor*etaz; - val_Proj_Jac_Tensor_i[0][2] = -factor*etay; - - val_Proj_Jac_Tensor_i[1][0] = -factor*etaz; - val_Proj_Jac_Tensor_i[1][1] = -factor*thetay; - val_Proj_Jac_Tensor_i[1][2] = -factor*etax; - - val_Proj_Jac_Tensor_i[2][0] = -factor*etay; - val_Proj_Jac_Tensor_i[2][1] = -factor*etax; - val_Proj_Jac_Tensor_i[2][2] = -factor*thetaz; - - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - val_Proj_Jac_Tensor_j[iVar][jVar] = -val_Proj_Jac_Tensor_i[iVar][jVar]; - - } - - if (nDim == 2) { - su2double thetax = theta + val_normal[0]*val_normal[0]/3.0; - su2double thetay = theta + val_normal[1]*val_normal[1]/3.0; - su2double etaz = val_normal[0]*val_normal[1]/3.0; - - - val_Proj_Jac_Tensor_i[0][0] = -factor*thetax; - val_Proj_Jac_Tensor_i[0][1] = -factor*etaz; - - val_Proj_Jac_Tensor_i[1][0] = -factor*etaz; - val_Proj_Jac_Tensor_i[1][1] = -factor*thetay; - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - val_Proj_Jac_Tensor_j[iVar][jVar] = -val_Proj_Jac_Tensor_i[iVar][jVar]; - } - -} diff --git a/SU2_CFD/src/numerics/pbflow.cpp b/SU2_CFD/src/numerics/pbflow.cpp deleted file mode 100644 index 7a940bea026..00000000000 --- a/SU2_CFD/src/numerics/pbflow.cpp +++ /dev/null @@ -1,573 +0,0 @@ -/*! - * \file pbflow.cpp - * \brief This file contains the numerical methods for incompressible flow. - * \author F. Palacios, T. Economon - * \version 6.0.1 "Falcon" - * - * The current SU2 release has been coordinated by the - * SU2 International Developers Society - * with selected contributions from the open-source community. - * - * The main research teams contributing to the current release are: - * - Prof. Juan J. Alonso's group at Stanford University. - * - Prof. Piero Colonna's group at Delft University of Technology. - * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. - * - Prof. Alberto Guardone's group at Polytechnic University of Milan. - * - Prof. Rafael Palacios' group at Imperial College London. - * - Prof. Vincent Terrapon's group at the University of Liege. - * - Prof. Edwin van der Weide's group at the University of Twente. - * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. - * - * Copyright 2012-2018, Francisco D. Palacios, Thomas D. Economon, - * Tim Albring, and the SU2 contributors. - * - * 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/pbflow.hpp" - -CUpwPB_Flow::CUpwPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - gravity = config->GetGravityForce(); - Froude = config->GetFroude(); - dynamic_grid = config->GetDynamic_Grid(); - - Diff_U = new su2double [nVar]; - Velocity_i = new su2double [nDim]; - Velocity_j = new su2double [nDim]; - Velocity_upw = new su2double [nDim]; - MeanVelocity = new su2double [nDim]; - Flux = new su2double [nDim]; - ProjFlux_i = new su2double [nVar]; - ProjFlux_j = new su2double [nVar]; - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - Jacobian_upw = new su2double* [nVar]; - - - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - Jacobian_upw[iVar] = new su2double [nVar]; - } - -} - -CUpwPB_Flow::~CUpwPB_Flow(void) { - - delete [] Diff_U; - delete [] Velocity_i; - delete [] Velocity_j; - delete [] Velocity_upw; - delete [] MeanVelocity; - delete [] Flux; - delete [] ProjFlux_i; - delete [] ProjFlux_j; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - delete [] Jacobian_upw[iVar]; - } - delete [] Jacobian_i; - delete [] Jacobian_j; - delete [] Jacobian_upw; - -} - -CNumerics::ResidualType<> CUpwPB_Flow::ComputeResidual(const CConfig *config) { - - - su2double MeanDensity, Flux0, Flux1, MeanPressure, Area, FF, Vel0, Vel1, ProjGridVelFlux = 0.0; - - - /*--- Primitive variables at point i and j ---*/ - Pressure_i = V_i[0]; Pressure_j = V_j[0]; - DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; - MeanDensity = 0.5*(DensityInc_i + DensityInc_j); - MeanPressure = 0.5*(Pressure_i + Pressure_j); - - Area = 0.0; - for(iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - - /*--- (rho*u_i) ---*/ - Face_Flux = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - Velocity_j[iDim] = V_j[iDim+1]; - MeanVelocity[iDim] = 0.5*(Velocity_i[iDim] + Velocity_j[iDim]); - Face_Flux += MeanDensity*MeanVelocity[iDim]*Normal[iDim]; - } - - if (dynamic_grid) { - ProjGridVelFlux = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjGridVelFlux += 0.5*MeanDensity*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - } - Face_Flux -= ProjGridVelFlux; - } - - /*--- Find upwind direction. ---*/ - Flux0 = 0.5*(Face_Flux + fabs(Face_Flux)) ; - Flux1 = 0.5*(Face_Flux - fabs(Face_Flux)) ; - - Upw_i = round(fabs(Flux0/(fabs(Face_Flux)+EPS))); - Upw_j = round(fabs(Flux1/(fabs(Face_Flux)+EPS))); - - /*--- Find flux. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Flux[iVar] = Flux0*V_i[iVar+1] + Flux1*V_j[iVar+1]; - Velocity_upw[iVar] = Upw_i*V_i[iVar+1] + Upw_j*V_j[iVar+1]; - if (dynamic_grid) Velocity_upw[iVar] -= (Upw_i*GridVel_i[iVar] + Upw_j*GridVel_j[iVar]); - } - - if (implicit) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_j[iVar][jVar] = 0.0; - Jacobian_i[iVar][jVar] = 0.0; - Jacobian_upw[iVar][jVar] = 0.0; - } - - GetInviscidPBProjJac(MeanDensity, Velocity_upw, Normal, 0.5, Jacobian_upw); - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = Upw_i*Jacobian_upw[iVar][jVar]; - Jacobian_j[iVar][jVar] = Upw_j*Jacobian_upw[iVar][jVar]; - } - } - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); -} - -CCentPB_Flow::CCentPB_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - gravity = config->GetGravityForce(); - Froude = config->GetFroude(); - dynamic_grid = config->GetDynamic_Grid(); - - Diff_U = new su2double [nVar]; - Velocity_i = new su2double [nDim]; - Velocity_j = new su2double [nDim]; - Velocity_upw = new su2double [nDim]; - MeanVelocity = new su2double [nDim]; - Flux = new su2double [nDim]; - ProjFlux_i = new su2double [nVar]; - ProjFlux_j = new su2double [nVar]; - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - Jacobian_upw = new su2double* [nVar]; - - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - Jacobian_upw[iVar] = new su2double [nVar]; - } -} - -CCentPB_Flow::~CCentPB_Flow(void) { - - delete [] Diff_U; - delete [] Velocity_i; - delete [] Velocity_j; - delete [] Velocity_upw; - delete [] MeanVelocity; - delete [] Flux; - delete [] ProjFlux_i; - delete [] ProjFlux_j; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - delete [] Jacobian_upw[iVar]; - } - delete [] Jacobian_i; - delete [] Jacobian_j; - delete [] Jacobian_upw; - -} - -CNumerics::ResidualType<> CCentPB_Flow::ComputeResidual(const CConfig *config) { - /*--- To do list - * 1. Find upwind direction, upwind node and downwind node - * 2. Compute normalised variable for the upwind node using the gradient of the upwind node - * 3. Find face velocity using central difference scheme - * 4. Use the ULTIMATE limiter to find the adjusted face velocity - * 5. Find residual as FaceFlux * AdjustedFaceVelocity - * */ - - - su2double MeanDensity, Flux0, Flux1, MeanPressure, Area, FF, Vel0, Vel1, ProjGridVelFlux = 0.0; - su2double dissipation, kappa=0.15, ProjVel_i = 0.0, ProjVel_j = 0.0; - - /*--- Conservative variables at point i and j ---*/ - - Pressure_i = V_i[0]; Pressure_j = V_j[0]; - DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; - - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - Velocity_j[iDim] = V_j[iDim+1]; - } - - Area = 0.0; - for(iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - - Face_Flux = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - Velocity_j[iDim] = V_j[iDim+1]; - MeanVelocity[iDim] = 0.5*(Velocity_i[iDim] + Velocity_j[iDim]); - Face_Flux += MeanDensity*MeanVelocity[iDim]*Normal[iDim]; - ProjVel_i += Velocity_i[iDim]*Normal[iDim]; - ProjVel_j += Velocity_j[iDim]*Normal[iDim]; - } - - if (dynamic_grid) { - ProjGridVelFlux = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjGridVelFlux += 0.5*MeanDensity*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - } - Face_Flux -= ProjGridVelFlux; - } - - Flux0 = 0.5*(Face_Flux + fabs(Face_Flux)) ; - Flux1 = 0.5*(Face_Flux - fabs(Face_Flux)) ; - - Upw_i = round(fabs(Flux0/(fabs(Face_Flux)+EPS))); - Upw_j = round(fabs(Flux1/(fabs(Face_Flux)+EPS))); - - for (iVar = 0; iVar < nVar; iVar++) { - Velocity_upw[iVar] = Upw_i*V_i[iVar+1] + Upw_j*V_j[iVar+1]; - if (dynamic_grid) Velocity_upw[iVar] -= (Upw_i*GridVel_i[iVar] + Upw_j*GridVel_j[iVar]); - Flux[iVar] = Face_Flux*MeanVelocity[iVar]; - } - - //Dissipation - su2double lambda_i = 0.0, lambda_j = 0.0; - lambda_i = 2.0*abs(ProjVel_i); - lambda_j = 2.0*abs(ProjVel_j); - - su2double lambda_mean = 0.5*(lambda_i + lambda_j); - if (lambda_mean < EPS) { - lambda_i = 2.0*abs(config->GetVelocity_Ref()*Area); - lambda_j = 2.0*abs(config->GetVelocity_Ref()*Area); - lambda_mean = abs(config->GetVelocity_Ref()*Area); - } - - su2double Param_p = 0.3, SF=0.0; - su2double Phi_i = pow((lambda_i/(4.0*lambda_mean)),Param_p); - su2double Phi_j = pow((lambda_j/(4.0*lambda_mean)),Param_p); - - if ((Phi_i + Phi_j) != 0.0) - SF = 4.0*Phi_i*Phi_j/(Phi_i + Phi_j); - - su2double sc0 = 3.0*(Neighbor_i + Neighbor_j)/(Neighbor_i*Neighbor_j); - su2double E_0 = kappa*sc0*2.0/3.0; - - su2double diss[MAXNDIM]; - for (iDim = 0; iDim < nDim; iDim++) - diss[iDim] = E_0*(DensityInc_i*Velocity_i[iDim] - DensityInc_j*Velocity_j[iDim])*SF*lambda_mean; - - for (iVar = 0; iVar < nVar; iVar++) - Flux[iVar] += diss[iVar]; - - /*--- For implicit schemes, compute jacobians based on the upwind scheme for stability issues. (See ANSYS user guide) ---*/ - if (implicit) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_j[iVar][jVar] = 0.0; - Jacobian_i[iVar][jVar] = 0.0; - Jacobian_upw[iVar][jVar] = 0.0; - } - - GetInviscidPBProjJac(MeanDensity, MeanVelocity, Normal, 0.5, Jacobian_upw); - //GetInviscidPBProjJac(&DensityInc_i, Velocity_upw, Normal, 0.5, val_Jacobian_upw); - /*GetInviscidPBProjJac(&DensityInc_i, Velocity_i, Normal, 0.5, Jacobian_i); - GetInviscidPBProjJac(&DensityInc_i, Velocity_j, Normal, 0.5, Jacobian_j);*/ - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = Upw_i*Jacobian_upw[iVar][jVar]; - Jacobian_j[iVar][jVar] = Upw_j*Jacobian_upw[iVar][jVar]; - } - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar][iVar] += E_0*DensityInc_i*SF*lambda_mean; - Jacobian_i[iVar][iVar] -= E_0*DensityInc_j*SF*lambda_mean; - } - - } - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); - -} - -CAvgGradPBInc_Flow::CAvgGradPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - - 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]; - } - - /*--- Primitive flow variables. ---*/ - - PrimVar_i = new su2double [nDim+4]; - PrimVar_j = new su2double [nDim+4]; - Mean_PrimVar = new su2double [nDim+4]; - - /*--- Incompressible flow, primitive variables nDim+2, (P, vx, vy, vz, rho) ---*/ - - Mean_GradPrimVar = new su2double*[nVar]; - - /*--- Incompressible flow, gradient primitive variables nDim+2, (P, vx, vy, vz, rho) ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - Mean_GradPrimVar[iVar] = new su2double[nDim]; - -} - -CAvgGradPBInc_Flow::~CAvgGradPBInc_Flow(void) { - - delete [] PrimVar_i; - delete [] PrimVar_j; - delete [] Mean_PrimVar; - - for (iVar = 0; iVar < nVar; iVar++) - delete [] Mean_GradPrimVar[iVar]; - delete [] Mean_GradPrimVar; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - } - delete [] Jacobian_i; - delete [] Jacobian_j; - -} - -CNumerics::ResidualType<> CAvgGradPBInc_Flow::ComputeResidual(const CConfig *config) { - - /*--- Normalized normal vector ---*/ - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; - - for (iVar = 0; iVar < nDim+4; iVar++) { - PrimVar_i[iVar] = V_i[iVar]; - PrimVar_j[iVar] = V_j[iVar]; - Mean_PrimVar[iVar] = 0.5*(PrimVar_i[iVar]+PrimVar_j[iVar]); - } - - /*--- Density and transport properties ---*/ - - Laminar_Viscosity_i = V_i[nDim+2]; Laminar_Viscosity_j = V_j[nDim+2]; - Eddy_Viscosity_i = V_i[nDim+3]; Eddy_Viscosity_j = V_j[nDim+3]; - - /*--- Mean transport properties ---*/ - - Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); - Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); - Mean_turb_ke = 0.0;//0.5*(turb_ke_i + turb_ke_j); - - /*--- Mean gradient approximation ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Mean_GradPrimVar[iVar][iDim] = 0.5*(PrimVar_Grad_i[iVar+1][iDim] + PrimVar_Grad_j[iVar+1][iDim]); - - /*--- Get projected flux tensor ---*/ - - GetViscousPBIncProjFlux(Mean_PrimVar, Mean_GradPrimVar, Normal, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, Mean_turb_ke); - - - /*--- Implicit part ---*/ - - if (implicit) { - - dist_ij = 0.0; proj_vector_ij = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - proj_vector_ij += (Coord_j[iDim]-Coord_i[iDim])*Normal[iDim]; - } - proj_vector_ij = proj_vector_ij/dist_ij; - dist_ij = sqrt(dist_ij); - - if (dist_ij == 0.0) { - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = 0.0; - Jacobian_j[iVar][jVar] = 0.0; - } - } - } - else { - GetViscousPBIncProjJacs(Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, dist_ij, UnitNormal, - Area, Jacobian_i, Jacobian_j); - } - - } - - return ResidualType<>(Proj_Flux_Tensor, Jacobian_i, Jacobian_j); -} - -CAvgGradCorrectedPBInc_Flow::CAvgGradCorrectedPBInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - - 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]; - } - - /*--- Primitive flow variables. ---*/ - - Mean_PrimVar = new su2double [nDim+4]; - PrimVar_i = new su2double [nDim+4]; - PrimVar_j = new su2double [nDim+4]; - Proj_Mean_GradPrimVar_Edge = new su2double [nVar]; - Edge_Vector = new su2double [nDim]; - - Mean_GradPrimVar = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Mean_GradPrimVar[iVar] = new su2double [nDim]; - -} - -CAvgGradCorrectedPBInc_Flow::~CAvgGradCorrectedPBInc_Flow(void) { - - delete [] Mean_PrimVar; - delete [] PrimVar_i; - delete [] PrimVar_j; - delete [] Proj_Mean_GradPrimVar_Edge; - delete [] Edge_Vector; - - for (iVar = 0; iVar < nVar; iVar++) - delete [] Mean_GradPrimVar[iVar]; - delete [] Mean_GradPrimVar; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - } - delete [] Jacobian_i; - delete [] Jacobian_j; - -} - -CNumerics::ResidualType<> CAvgGradCorrectedPBInc_Flow::ComputeResidual(const CConfig *config) { - - unsigned short nPrimVarGrad = nDim+2, nPrimVar = nDim+4; - - /*--- Normalized normal vector ---*/ - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; - - /*--- Conversion to Primitive Variables (P, u, v, w, rho, mu, muT) ---*/ - - for (iVar = 0; iVar < nPrimVar; iVar++) { - PrimVar_i[iVar] = V_i[iVar]; - PrimVar_j[iVar] = V_j[iVar]; - Mean_PrimVar[iVar] = 0.5*(PrimVar_i[iVar]+PrimVar_j[iVar]); - } - - /*--- Density and transport properties ---*/ - - DensityInc_i = V_i[nDim+1]; DensityInc_j = V_j[nDim+1]; - Laminar_Viscosity_i = V_i[nDim+2]; Laminar_Viscosity_j = V_j[nDim+2]; - Eddy_Viscosity_i = V_i[nDim+3]; Eddy_Viscosity_j = V_j[nDim+3]; - - /*--- Mean transport properties ---*/ - - Mean_Laminar_Viscosity = 0.5*(Laminar_Viscosity_i + Laminar_Viscosity_j); - Mean_Eddy_Viscosity = 0.5*(Eddy_Viscosity_i + Eddy_Viscosity_j); - Mean_turb_ke = 0.0;//0.5*(turb_ke_i + turb_ke_j); - - /*--- Compute vector going from iPoint to jPoint ---*/ - - dist_ij_2 = 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]; - } - - /*--- Projection of the mean gradient in the direction of the edge ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradPrimVar_Edge[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPrimVar[iVar][iDim] = 0.5*(PrimVar_Grad_i[iVar+1][iDim] + PrimVar_Grad_j[iVar+1][iDim]); - Proj_Mean_GradPrimVar_Edge[iVar] += Mean_GradPrimVar[iVar][iDim]*Edge_Vector[iDim]; - } - if (dist_ij_2 != 0.0) { - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPrimVar[iVar][iDim] -= (Proj_Mean_GradPrimVar_Edge[iVar] - - (PrimVar_j[iVar+1]-PrimVar_i[iVar+1]))*Edge_Vector[iDim] / dist_ij_2; - } - } - } - - /*--- Get projected flux tensor ---*/ - - GetViscousPBIncProjFlux(Mean_PrimVar, Mean_GradPrimVar, Normal, Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, Mean_turb_ke); - - /*--- Implicit part ---*/ - - if (implicit) { - - proj_vector_ij = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - proj_vector_ij += (Coord_j[iDim]-Coord_i[iDim])*Normal[iDim]; - } - proj_vector_ij = proj_vector_ij/dist_ij_2; - - if (dist_ij_2 == 0.0) { - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = 0.0; - Jacobian_j[iVar][jVar] = 0.0; - } - } - } - else { - GetViscousPBIncProjJacs(Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, sqrt(dist_ij_2), UnitNormal, - Area, Jacobian_i, Jacobian_j); - } - - } - - return ResidualType<>(Proj_Flux_Tensor, Jacobian_i, Jacobian_j); -} diff --git a/SU2_CFD/src/numerics/poisson.cpp b/SU2_CFD/src/numerics/poisson.cpp deleted file mode 100644 index f096a1b70fc..00000000000 --- a/SU2_CFD/src/numerics/poisson.cpp +++ /dev/null @@ -1,421 +0,0 @@ -/*! - * \file numerics_direct_poisson.cpp - * \brief This file contains all the convective term discretization. - * \author F. Palacios - * \version 6.1.0 "Falcon" - * - * The current SU2 release has been coordinated by the - * SU2 International Developers Society - * with selected contributions from the open-source community. - * - * The main research teams contributing to the current release are: - * - Prof. Juan J. Alonso's group at Stanford University. - * - Prof. Piero Colonna's group at Delft University of Technology. - * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. - * - Prof. Alberto Guardone's group at Polytechnic University of Milan. - * - Prof. Rafael Palacios' group at Imperial College London. - * - Prof. Vincent Terrapon's group at the University of Liege. - * - Prof. Edwin van der Weide's group at the University of Twente. - * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. - * - * Copyright 2012-2018, Francisco D. Palacios, Thomas D. Economon, - * Tim Albring, and the SU2 contributors. - * - * 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/poisson.hpp" - - -CAvgGradCorrected_Poisson::CAvgGradCorrected_Poisson(unsigned short val_nDim, unsigned short val_nVar, - CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = false; - implicit = (config->GetKind_TimeIntScheme_Poisson() == EULER_IMPLICIT); - direct = false; - - - Edge_Vector = new su2double [nDim]; - Proj_Mean_GradPoissonVar_Edge = new su2double [nVar]; - Proj_Mean_GradPoissonVar_Kappa = new su2double [nVar]; - Proj_Mean_GradPoissonVar_Corrected = new su2double [nVar]; - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - Mean_GradPoissonVar = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nDim]; - Jacobian_j[iVar] = new su2double [nDim]; - Mean_GradPoissonVar[iVar] = new su2double [nDim]; - } - - Poisson_Coeff_i = 1.0; - Poisson_Coeff_j = 1.0; - - Mom_Coeff_i = new su2double [nDim]; - Mom_Coeff_j = new su2double [nDim]; - -} - - -CAvgGradCorrected_Poisson::~CAvgGradCorrected_Poisson(void) { - - delete [] Edge_Vector; - delete [] Proj_Mean_GradPoissonVar_Edge; - delete [] Proj_Mean_GradPoissonVar_Kappa; - delete [] Proj_Mean_GradPoissonVar_Corrected; - delete [] Mom_Coeff_i; - delete [] Mom_Coeff_j; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Mean_GradPoissonVar[iVar]; - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - } - delete [] Mean_GradPoissonVar; - delete [] Jacobian_i; - delete [] Jacobian_j; -} - -CNumerics::ResidualType<> CAvgGradCorrected_Poisson::ComputeResidual(const CConfig *config) { - - su2double Coeff_Mean; - su2double *MomCoeffxNormal = new su2double[nDim]; - su2double *MomCoeffxNormalCorrected = new su2double[nDim]; - su2double Mean_GradPoissonVar_Edge[MAXNDIM], GradPoisson[MAXNDIM]; - su2double Mean_GradPoissonVar_Face[MAXNDIM][MAXNDIM], Num, Den; - - Poisson_Coeff_Mean = 1.0;//0.5*(Poisson_Coeff_i + Poisson_Coeff_j); - - /*--- Multiply the normal with the coefficients of the momentum equation - * and use it instead of the normal during the projection operation ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - Coeff_Mean = 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]) ; - MomCoeffxNormal[iDim] = Coeff_Mean*Normal[iDim]; - } - - /*--- Compute vector going from iPoint to jPoint ---*/ - /*--- Minimum correction approach. ---*/ - /*dist_ij_2 = 0; proj_vector_ij = 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]; - proj_vector_ij += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2;*/ - - /*--- Over relaxed approach. ---*/ - dist_ij_2 = 0; proj_vector_ij = 0; - Num = 0.0; Den = 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]; - Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; - Den += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; - } - if (Den == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = Num/Den; - - /*--- Orthogonal correction approach. ---*/ - /*dist_ij_2 = 0; proj_vector_ij = 0; - Num = 0.0; Den = 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]; - Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; - } - Den = dist_ij_2; - if (Den == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = sqrt(Num/Den);*/ - - /*--- Correct the face gradient for odd-even decoupling ---*/ - /*--- Steps are - * 1. Interpolate the gradient at the face -> du/ds|_in - * 2. Find the projection of the interpolated gradient on the edge vector -> (du/ds|_in) . e - * 3. Find the gradient as the difference between the neighboring nodes -> (u_j - u_i)/ds - * 4. Correct the gradient at the face using du/ds = du/ds|_in + ( (u_j - u_i)/ds - (du/ds|_in . e) ) . e ---*/ - - /*for (iVar = 0; iVar < nVar; iVar++) { - Mean_GradPoissonVar_Edge[iVar] = 0.0; - - GradPoisson[iVar] = (Poissonval_j - Poissonval_i)/sqrt(dist_ij_2); - - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); - - Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]/dist_ij_2; - } - } - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPoissonVar_Face[iVar][iDim] = Mean_GradPoissonVar[iVar][iDim] + - (GradPoisson[iVar] - Mean_GradPoissonVar_Edge[iVar])*Edge_Vector[iDim]/dist_ij_2; - } - }*/ - - - /*--- Mean gradient approximation. Projection of the mean gradient - in the direction of the edge ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradPoissonVar_Edge[iVar] = 0.0; - Proj_Mean_GradPoissonVar_Kappa[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); - Proj_Mean_GradPoissonVar_Kappa[iVar] += Mean_GradPoissonVar[iVar][iDim]*MomCoeffxNormal[iDim]; - Proj_Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]; - } - Proj_Mean_GradPoissonVar_Corrected[iVar] = Proj_Mean_GradPoissonVar_Kappa[iVar]; - Proj_Mean_GradPoissonVar_Corrected[iVar] -= Proj_Mean_GradPoissonVar_Edge[iVar]*proj_vector_ij - - (Poissonval_j-Poissonval_i)*proj_vector_ij; - } - - /*--- Jacobians for implicit scheme ---*/ - - if (implicit) { - Jacobian_i[0][0] = -Poisson_Coeff_Mean*proj_vector_ij; - Jacobian_j[0][0] = Poisson_Coeff_Mean*proj_vector_ij; - } - - /*AD::SetPreaccOut(val_residual, nVar); - AD::EndPreacc();*/ - - delete [] MomCoeffxNormal; - delete [] MomCoeffxNormalCorrected; - - return ResidualType<>(Proj_Mean_GradPoissonVar_Corrected, Jacobian_i, Jacobian_j); -} - -CAvgGrad_Poisson::CAvgGrad_Poisson(unsigned short val_nDim, unsigned short val_nVar, - CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = false; - implicit = (config->GetKind_TimeIntScheme_Poisson() == EULER_IMPLICIT); // TODO: PBFlow: FIXED - direct = false; - - - Edge_Vector = new su2double [nDim]; - Proj_Mean_GradPoissonVar_Normal = new su2double [nVar]; - Proj_Mean_GradPoissonVar_Corrected = new su2double [nVar]; - Mom_Coeff_i = new su2double [nDim]; - Mom_Coeff_j = new su2double [nDim]; - - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - Mean_GradPoissonVar = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nDim]; - Jacobian_j[iVar] = new su2double [nDim]; - Mean_GradPoissonVar[iVar] = new su2double [nDim]; - } - - Poisson_Coeff_i = 1.0; - Poisson_Coeff_j = 1.0; -} - - -CAvgGrad_Poisson::~CAvgGrad_Poisson(void) { - - unsigned int iDim; - delete [] Edge_Vector; - delete [] Proj_Mean_GradPoissonVar_Normal; - delete [] Proj_Mean_GradPoissonVar_Corrected; - delete [] Mom_Coeff_i; - delete [] Mom_Coeff_j; - - for (iVar = 0; iVar < nVar; iVar++) { - delete [] Mean_GradPoissonVar[iVar]; - delete [] Jacobian_i[iVar]; - delete [] Jacobian_j[iVar]; - } - delete [] Mean_GradPoissonVar; - delete [] Jacobian_i; - delete [] Jacobian_j; - -} - -CNumerics::ResidualType<> CAvgGrad_Poisson::ComputeResidual(const CConfig *config) { - - - su2double Coeff_Mean, Num, Den; - su2double *MomCoeffxNormal = new su2double[nDim]; - su2double Mean_GradPoissonVar_Edge[MAXNDIM], GradPoisson[MAXNDIM]; - su2double Mean_GradPoissonVar_Face[MAXNDIM][MAXNDIM]; - - - - Poisson_Coeff_Mean = 1.0;//0.5*(Poisson_Coeff_i + Poisson_Coeff_j); - - /*--- Multiply the normal with the coefficients of the momentum equation - * and use it instead of the normal during the projection operation ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - Coeff_Mean = 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]) ; - MomCoeffxNormal[iDim] = Coeff_Mean*Normal[iDim]; - } - - /*--- Compute vector going from iPoint to jPoint ---*/ - /*--- Minimum correction approach. ---*/ - /*dist_ij_2 = 0; proj_vector_ij = 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]; - proj_vector_ij += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2;*/ - - /*--- Over relaxed approach. ---*/ - dist_ij_2 = 0; proj_vector_ij = 0; - Num = 0.0; Den = 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]; - Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; - Den += Edge_Vector[iDim]*MomCoeffxNormal[iDim]; - } - if (Den == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = Num/Den; - - /*--- Orthogonal correction approach. ---*/ - /*dist_ij_2 = 0; proj_vector_ij = 0; - Num = 0.0; Den = 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]; - Num += MomCoeffxNormal[iDim]*MomCoeffxNormal[iDim]; - } - Den = dist_ij_2; - if (Den == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = sqrt(Num/Den);*/ - - /*--- Correct the face gradient for odd-even decoupling ---*/ - /*--- Steps are - * 1. Interpolate the gradient at the face -> du/ds|_in - * 2. Find the projection of the interpolated gradient on the edge vector -> (du/ds|_in) . e - * 3. Find the gradient as the difference between the neighboring nodes -> (u_j - u_i)/ds - * 4. Correct the gradient at the face using du/ds = du/ds|_in + ( (u_j - u_i)/ds - (du/ds|_in . e) ) . e ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Mean_GradPoissonVar_Edge[iVar] = 0.0; - - GradPoisson[iVar] = (Poissonval_j - Poissonval_i)/sqrt(dist_ij_2); - - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPoissonVar[iVar][iDim] = 0.5*(ConsVar_Grad_i[iVar][iDim] + ConsVar_Grad_j[iVar][iDim]); - - Mean_GradPoissonVar_Edge[iVar] += Mean_GradPoissonVar[iVar][iDim]*Edge_Vector[iDim]/sqrt(dist_ij_2); - } - } - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Mean_GradPoissonVar_Face[iVar][iDim] = Mean_GradPoissonVar[iVar][iDim] + - (GradPoisson[iVar] - Mean_GradPoissonVar_Edge[iVar])*Edge_Vector[iDim]/sqrt(dist_ij_2); - } - } - - /*--- Mean gradient approximation. Projection of the mean gradient in the direction of the edge ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Proj_Mean_GradPoissonVar_Normal[iVar] = 0.0; - Proj_Mean_GradPoissonVar_Corrected[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Proj_Mean_GradPoissonVar_Normal[iVar] += Mean_GradPoissonVar_Face[iVar][iDim]*MomCoeffxNormal[iDim]; - } - } - - /*--- Jacobians for implicit scheme ---*/ - if (implicit) { - Jacobian_i[0][0] = -Poisson_Coeff_Mean*proj_vector_ij; - Jacobian_j[0][0] = Poisson_Coeff_Mean*proj_vector_ij; - } - - delete [] MomCoeffxNormal; - - return ResidualType<>(Proj_Mean_GradPoissonVar_Normal, Jacobian_i, Jacobian_j); -} - - - -CPressure_Poisson::CPressure_Poisson(unsigned short val_nDim, unsigned short val_nVar, - CConfig *config) : CNumerics(val_nDim, val_nVar, config) { - - implicit = true; - - Edge_Vector = new su2double [nDim]; - - Poisson_Coeff_i = 1.0; - Poisson_Coeff_j = 1.0; - - Mom_Coeff_i = new su2double [nDim]; - Mom_Coeff_j = new su2double [nDim]; - -} - - -CPressure_Poisson::~CPressure_Poisson(void) { - - unsigned int iDim; - delete [] Edge_Vector; - - delete Mom_Coeff_i; - delete Mom_Coeff_j; -} - -void CPressure_Poisson::ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) { - - - su2double Area; - su2double UnitNormal[MAXNDIM]; - - /*--- Compute vector going from iPoint to jPoint ---*/ - - dist_ij_2 = 0; proj_vector_ij = 0; - Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - UnitNormal[iDim] = Normal[iDim]/Area; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; - - - Poisson_Coeff_Mean = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Poisson_Coeff_Mean += 0.5*(Mom_Coeff_i[iDim] + Mom_Coeff_j[iDim]); - - - val_residual[0] = 0.5*Poisson_Coeff_Mean*(Poissonval_i + Poissonval_j); - - if (implicit){ - Jacobian_i[0][0] = Poisson_Coeff_Mean; - Jacobian_j[0][0] = -Poisson_Coeff_Mean; - } -} - - -CSource_PoissonFVM::CSource_PoissonFVM(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { } - -CSource_PoissonFVM::~CSource_PoissonFVM(void) { } - -void CSource_PoissonFVM::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { - - if (config->GetKind_Incomp_System()==ENUM_INCOMP_SYSTEM::PRESSURE_BASED) - val_residual[0] = Source_Term; - else - val_residual[0] = Source_Term*Volume; -} From d973f1997ea8ed47a70eda22a2dad498ad7e8e4a Mon Sep 17 00:00:00 2001 From: Thijs Aalbers Date: Tue, 12 May 2026 13:46:17 +0200 Subject: [PATCH 4/4] Introduce pressure-based solver architecture in minimal form Implements the full class structure required for the pressure-based solver in a minimal form. The code compiles and runs but does not yet contain any numerical/physical implementation. Future work will focus on implementing solver logic. Note: In the previous attempts (see related work) there is noticeable code duplication between CIncEuler and CPBIncEuler. The final architecture may be revised depending on how the implementation of CPBIncEuler evolves. --- .../include/iteration/CPBFluidIteration.hpp | 75 ++++++ SU2_CFD/include/solvers/CPBIncEulerSolver.hpp | 72 ++++++ SU2_CFD/include/solvers/CPBIncNSSolver.hpp | 47 ++++ SU2_CFD/include/solvers/CPoissonSolver.hpp | 119 +++++++++ SU2_CFD/include/solvers/CSolverFactory.hpp | 1 + .../include/variables/CPBIncEulerVariable.hpp | 80 ++++++ .../include/variables/CPBIncNSVariable.hpp | 59 +++++ .../include/variables/CPoissonVariable.hpp | 59 +++++ SU2_CFD/src/drivers/CDriver.cpp | 9 +- SU2_CFD/src/iteration/CIterationFactory.cpp | 6 + SU2_CFD/src/iteration/CPBFluidIteration.cpp | 64 +++++ SU2_CFD/src/meson.build | 7 + SU2_CFD/src/solvers/CPBIncEulerSolver.cpp | 232 ++++++++++++++++++ SU2_CFD/src/solvers/CPBIncNSSolver.cpp | 42 ++++ SU2_CFD/src/solvers/CPoissonSolver.cpp | 41 ++++ SU2_CFD/src/solvers/CSolverFactory.cpp | 38 ++- SU2_CFD/src/variables/CPBIncEulerVariable.cpp | 37 +++ SU2_CFD/src/variables/CPBIncNSVariable.cpp | 38 +++ SU2_CFD/src/variables/CPoissonVariable.cpp | 0 19 files changed, 1020 insertions(+), 6 deletions(-) create mode 100644 SU2_CFD/include/iteration/CPBFluidIteration.hpp create mode 100644 SU2_CFD/include/solvers/CPBIncEulerSolver.hpp create mode 100644 SU2_CFD/include/solvers/CPBIncNSSolver.hpp create mode 100644 SU2_CFD/include/solvers/CPoissonSolver.hpp create mode 100644 SU2_CFD/include/variables/CPBIncEulerVariable.hpp create mode 100644 SU2_CFD/include/variables/CPBIncNSVariable.hpp create mode 100644 SU2_CFD/include/variables/CPoissonVariable.hpp create mode 100644 SU2_CFD/src/iteration/CPBFluidIteration.cpp create mode 100644 SU2_CFD/src/solvers/CPBIncEulerSolver.cpp create mode 100644 SU2_CFD/src/solvers/CPBIncNSSolver.cpp create mode 100644 SU2_CFD/src/solvers/CPoissonSolver.cpp create mode 100644 SU2_CFD/src/variables/CPBIncEulerVariable.cpp create mode 100644 SU2_CFD/src/variables/CPBIncNSVariable.cpp create mode 100644 SU2_CFD/src/variables/CPoissonVariable.cpp diff --git a/SU2_CFD/include/iteration/CPBFluidIteration.hpp b/SU2_CFD/include/iteration/CPBFluidIteration.hpp new file mode 100644 index 00000000000..a102b9fa2c8 --- /dev/null +++ b/SU2_CFD/include/iteration/CPBFluidIteration.hpp @@ -0,0 +1,75 @@ +/*! + * \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 Preprocessing to prepare for an iteration of the physics. + * \param[in] ??? - Description here. + */ + void Preprocess(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; + + /*! + * \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/solvers/CPBIncEulerSolver.hpp b/SU2_CFD/include/solvers/CPBIncEulerSolver.hpp new file mode 100644 index 00000000000..45add069e71 --- /dev/null +++ b/SU2_CFD/include/solvers/CPBIncEulerSolver.hpp @@ -0,0 +1,72 @@ +/*! + * \file CPBIncEulerSolver.hpp + * \brief Headers of the CPBIncEulerSolver class + * \author F. Palacios, T. Economon, T. Albring + * \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 "CFVMFlowSolverBase.hpp" +#include "../variables/CPBIncEulerVariable.hpp" + +class CPBIncEulerSolver : public CFVMFlowSolverBase { +protected: + +public: + /*! + * \brief Constructor of the class. + */ + //CPBIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) {}; + + CPBIncEulerSolver() = delete; + + /*! + * \brief Constructor of the class. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Grid level. + * \param[in] navier_stokes - True when the constructor is called by the derived class CPBIncNSSolver. + */ + CPBIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, const bool navier_stokes = false); + + /*! + * \brief Destructor of the class. + */ + ~CPBIncEulerSolver(void) override {} + + /*! + * \brief Set reference values for pressure, forces, etc. + */ + void SetReferenceValues(const CConfig& config) final {} + + /*! + * \brief Print verification error to screen. + * \param[in] config - Definition of the particular problem. + */ + void PrintVerificationError(const CConfig* config) const final {} + + + +}; diff --git a/SU2_CFD/include/solvers/CPBIncNSSolver.hpp b/SU2_CFD/include/solvers/CPBIncNSSolver.hpp new file mode 100644 index 00000000000..c20f69e626a --- /dev/null +++ b/SU2_CFD/include/solvers/CPBIncNSSolver.hpp @@ -0,0 +1,47 @@ +/*! + * \file CPBIncNSSolver.hpp + * \brief Headers of the CPBIncNSSolver class + * \author F. Palacios, T. Economon, T. Albring + * \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 "CPBIncEulerSolver.hpp" + +/*! + * \class CPBIncNSSolver + * \brief Main class for defining the pressure based incompressible Navier-Stokes flow solver. + * \ingroup Navier_Stokes_Equations + * \author F. Palacios, T. Economon, T. Albring + */ +class CPBIncNSSolver final : public CPBIncEulerSolver { +public: + + /*! + * \brief Constructor of the class. + */ + CPBIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh); + + +}; diff --git a/SU2_CFD/include/solvers/CPoissonSolver.hpp b/SU2_CFD/include/solvers/CPoissonSolver.hpp new file mode 100644 index 00000000000..49a9f723e8f --- /dev/null +++ b/SU2_CFD/include/solvers/CPoissonSolver.hpp @@ -0,0 +1,119 @@ +/*! + * \file CPoissonSolver.hpp + * \brief Headers of the CPoissonSolver 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 "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: + +public: + + /*! + * \brief Constructor of the class. + */ + // CPoissonSolver(void); + /* + * \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 Destructor of the class. + */ + // ~CPoissonSolver(void); + + /*! + * \brief Load a solution from a restart file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all of the solvers. + * \param[in] config - Definition of the particular problem. + * \param[in] val_iter - Current external iteration number. + * \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 {} + + /*! + * \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. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] Iteration - Index of the current iteration. + */ + void SetTime_Step(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned long Iteration) override {} + + /*! + * \brief Set the initial condition for the FEM structural problem. + * \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] ExtIter - External iteration. + */ + void SetInitialCondition(CGeometry **geometry, + CSolver ***solver_container, + CConfig *config, + unsigned long TimeIter) override {} + + /*! + * \brief Set the total residual adding the term that comes from the Dual Time-Stepping Strategy. + * \param[in] geometry - Geometric definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + */ + void SetResidual_DualTime(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iRKStep, + unsigned short iMesh, + unsigned short RunTime_EqSystem) override {} + + + + +}; 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/CPBIncEulerVariable.hpp b/SU2_CFD/include/variables/CPBIncEulerVariable.hpp new file mode 100644 index 00000000000..b0f098189b8 --- /dev/null +++ b/SU2_CFD/include/variables/CPBIncEulerVariable.hpp @@ -0,0 +1,80 @@ +/*! + * \file CPBIncEulerVariable.hpp + * \brief Class for defining the variables of the pressure based incompressible Euler 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 +#include "CFlowVariable.hpp" + +/*! + * \class CPBIncEulerVariable + * \brief Class for defining the variables of the pressure based incompressible Euler solver. + * \ingroup Euler_Equations + * \author F. Palacios, T. Economon, T. Albring + */ +class CPBIncEulerVariable : public CFlowVariable { +public: + static constexpr size_t MAXNVAR = 13; + + template + struct CIndices { + const IndexType nDim; + CIndices(IndexType ndim, IndexType) : nDim(ndim) {} + inline IndexType NDim() const { return nDim; } + inline IndexType NSpecies() const { return 0; } + inline IndexType Pressure() const { return 0; } + inline IndexType Velocity() const { return 1; } + inline IndexType Temperature() const { return nDim+1; } + inline IndexType Density() const { return nDim+2; } + inline IndexType Enthalpy() const { return nDim+3; } + inline IndexType Beta() const { return nDim+4; } + inline IndexType SoundSpeed() const { return Beta(); } + inline IndexType LaminarViscosity() const { return nDim+5; } + inline IndexType EddyViscosity() const { return nDim+6; } + inline IndexType ThermalConductivity() const { return nDim+7; } + inline IndexType CpTotal() const { return nDim+8; } + inline IndexType CvTotal() const { return nDim+9; } + + /*--- For compatible interface with NEMO. ---*/ + inline IndexType SpeciesDensities() const { return std::numeric_limits::max(); } + inline IndexType Temperature_ve() const { return std::numeric_limits::max(); } + }; + +protected: + const CIndices indices; + +public: + + CPBIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy, + unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config); + + inline CMatrixView GetVelocityGradient(unsigned long iPoint) const final { + return Gradient_Primitive(iPoint, indices.Velocity()); + } + + +}; diff --git a/SU2_CFD/include/variables/CPBIncNSVariable.hpp b/SU2_CFD/include/variables/CPBIncNSVariable.hpp new file mode 100644 index 00000000000..d9837bd9d29 --- /dev/null +++ b/SU2_CFD/include/variables/CPBIncNSVariable.hpp @@ -0,0 +1,59 @@ +/*! + * \file CPBIncNSVariable.hpp + * \brief Class for defining the variables of the pressure based incompressible + Navier-Stokes 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 "CPBIncEulerVariable.hpp" + +/*! + * \class CPBIncNSVariable + * \brief Class for defining the variables of the incompressible Navier-Stokes solver. + * \ingroup Navier_Stokes_Equations + * \author F. Palacios, T. Economon, T. Albring + */ +class CPBIncNSVariable final : public CPBIncEulerVariable { +private: + const bool Energy; /*!< \brief Flag for Energy equation in incompressible flows. */ + +public: + + /*! + * \brief Constructor of the class. + * \param[in] pressure - value of the pressure. + * \param[in] velocity - Value of the flow velocity (initialization value). + * \param[in] temperature - Value of the temperature (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. + */ + CPBIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature, + unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config); + + +}; diff --git a/SU2_CFD/include/variables/CPoissonVariable.hpp b/SU2_CFD/include/variables/CPoissonVariable.hpp new file mode 100644 index 00000000000..1d730ed5130 --- /dev/null +++ b/SU2_CFD/include/variables/CPoissonVariable.hpp @@ -0,0 +1,59 @@ +/*! + * \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 { +public: + static constexpr size_t MAXNVAR = 1; /*!< \brief Max number of variables, for static arrays. */ + + /*! + * \brief Constructor of the class. + * \param[in] heat - 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 heat, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) + : CScalarVariable(npoint, ndim, nvar, config) {} + + /*! + * \brief Get the temperature of the point. + * \return Value of the temperature of the point. + */ + // inline su2double GetTemperature(unsigned long iPoint) const final { return Solution(iPoint, 0); } + +}; diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index b2f4ec1360a..9667fdb9393 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -51,6 +51,7 @@ #include "../../include/variables/CEulerVariable.hpp" #include "../../include/variables/CIncEulerVariable.hpp" +#include "../../include/variables/CPBIncEulerVariable.hpp" #include "../../include/variables/CNEMOEulerVariable.hpp" #include "../../include/numerics/template.hpp" @@ -2083,10 +2084,10 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver if (incompressible) if (!pressure_based) InstantiateTurbulentNumerics >(nVar_Turb, offset, config, - solver[MESH_0][TURB_SOL], numerics); - // else - // InstantiateTurbulentNumerics >(nVar_Turb, offset, config, - // solver[MESH_0][TURB_SOL], numerics); + solver[MESH_0][TURB_SOL], numerics); + else + InstantiateTurbulentNumerics >(nVar_Turb, offset, config, + solver[MESH_0][TURB_SOL], numerics); else if (NEMO_ns) InstantiateTurbulentNumerics >(nVar_Turb, offset, config, diff --git a/SU2_CFD/src/iteration/CIterationFactory.cpp b/SU2_CFD/src/iteration/CIterationFactory.cpp index 761ceca98f8..5b5b084fa60 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() == ENUM_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; diff --git a/SU2_CFD/src/iteration/CPBFluidIteration.cpp b/SU2_CFD/src/iteration/CPBFluidIteration.cpp new file mode 100644 index 00000000000..9c806e5c5d7 --- /dev/null +++ b/SU2_CFD/src/iteration/CPBFluidIteration.cpp @@ -0,0 +1,64 @@ +/*! + * \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" + +void CPBFluidIteration::Preprocess(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) { } + +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 unsteady = (config[val_iZone]->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config[val_iZone]->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); + 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()); + + /*--- Setting up iteration values depending on if this is a + steady or an unsteady simulation */ + + const auto InnerIter = config[val_iZone]->GetInnerIter(); + const auto TimeIter = config[val_iZone]->GetTimeIter(); + + /*--- Solve the Euler, Navier-Stokes, RANS equations. ---*/ + + /*--- Solve the momentum equations. ---*/ + + /*--- Solve the pressure correction (poisson) equation ---*/ + + /*--- Correct pressure and velocities. ---*/ + + /*--- T.B.D. ---*/ + +} diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3b40822a34f..8aee51243dc 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -55,6 +55,7 @@ su2_cfd_src += files(['output/COutputFactory.cpp', 'output/tools/CWindowingTools.cpp']) su2_cfd_src += files(['variables/CIncNSVariable.cpp', + 'variables/CPBIncNSVariable.cpp', 'variables/CTransLMVariable.cpp', 'variables/CAdjEulerVariable.cpp', 'variables/CHeatVariable.cpp', @@ -81,6 +82,8 @@ su2_cfd_src += files(['variables/CIncNSVariable.cpp', 'variables/CAdjTurbVariable.cpp', 'variables/CFlowVariable.cpp', 'variables/CIncEulerVariable.cpp', + 'variables/CPBIncEulerVariable.cpp', + 'variables/CPoissonVariable.cpp', 'variables/CEulerVariable.cpp', 'variables/CNEMOEulerVariable.cpp', 'variables/CNEMONSVariable.cpp', @@ -103,7 +106,10 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CFEM_DG_NSSolver.cpp', 'solvers/CHeatSolver.cpp', 'solvers/CIncEulerSolver.cpp', + 'solvers/CPBIncEulerSolver.cpp', 'solvers/CIncNSSolver.cpp', + 'solvers/CPBIncNSSolver.cpp', + 'solvers/CPoissonSolver.cpp', 'solvers/CMeshSolver.cpp', 'solvers/CNEMOEulerSolver.cpp', 'solvers/CNEMONSSolver.cpp', @@ -184,6 +190,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/solvers/CPBIncEulerSolver.cpp b/SU2_CFD/src/solvers/CPBIncEulerSolver.cpp new file mode 100644 index 00000000000..b402713ceec --- /dev/null +++ b/SU2_CFD/src/solvers/CPBIncEulerSolver.cpp @@ -0,0 +1,232 @@ +/*! + * \file CPBIncEulerSolver.cpp + * \brief Main subroutines for solving pressure based incompressible flow (Euler, Navier-Stokes, etc.). + * \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/CPBIncEulerSolver.hpp" +#include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../include/fluid/CConstantDensity.hpp" +#include "../../include/fluid/CIncIdealGas.hpp" +#include "../../include/fluid/CIncIdealGasPolynomial.hpp" +#include "../../include/variables/CPBIncNSVariable.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../include/fluid/CFluidScalar.hpp" +#include "../../include/fluid/CFluidFlamelet.hpp" +#include "../../include/fluid/CFluidModel.hpp" + +// Temporarily directly copied from CIncEulerSolver's constructor to adhere to necessary +// initializations, +CPBIncEulerSolver::CPBIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, + const bool navier_stokes) : + CFVMFlowSolverBase(*geometry, *config) { + SU2_ZONE_SCOPED + + + /*--- Based on the navier_stokes boolean, determine if this constructor is + * being called by itself, or by its derived class CIncNSSolver. ---*/ + const string description = navier_stokes? "Navier-Stokes" : "Euler"; + + unsigned short iMarker; + ifstream restart_file; + bool restart = (config->GetRestart() || config->GetRestart_Flow()); + int Unst_RestartIter = 0; + bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)); + bool time_stepping = config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING; + bool adjoint = (config->GetContinuous_Adjoint()) || (config->GetDiscrete_Adjoint()); + const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED; + const su2double* scalar_init = nullptr; + + /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ + dynamic_grid = config->GetDynamic_Grid(); + + /*--- Store the multigrid level. ---*/ + MGLevel = iMesh; + + /*--- Check for a restart file to evaluate if there is a change in the angle of attack + before computing all the non-dimesional quantities. ---*/ + + if (restart) { + + if (iMesh == MESH_0) { + + /*--- Modify file name for a dual-time unsteady restart ---*/ + + if (dual_time) { + if (adjoint) Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; + else if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) + Unst_RestartIter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + else Unst_RestartIter = SU2_TYPE::Int(config->GetRestart_Iter())-2; + } + + /*--- Modify file name for a time stepping unsteady restart ---*/ + + if (time_stepping) { + if (adjoint) Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; + else Unst_RestartIter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + } + + } + + if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) { + if (rank==MASTER_NODE) cout << "Setting streamwise periodic pressure drop from restart metadata file." << endl; + } + + auto filename_ = config->GetFilename("flow", ".meta", Unst_RestartIter); + Read_SU2_Restart_Metadata(geometry, config, adjoint, filename_); + } + + + /*--- Set the gamma value ---*/ + + Gamma = config->GetGamma(); + Gamma_Minus_One = Gamma - 1.0; + + /*--- Define geometry constants in the solver structure. + * Incompressible flow, primitive variables (P, vx, vy, vz, T, rho, beta, lamMu, EddyMu, Kt_eff, Cp, Cv) ---*/ + + nDim = geometry->GetnDim(); + + /*--- Make sure to align the sizes with the constructor of CIncEulerVariable. ---*/ + nVar = nDim + 2; + nPrimVar = nDim + 10; + /*--- Centered schemes only need gradients for viscous fluxes (T and v, but we need also to include P). ---*/ + nPrimVarGrad = nDim + (centered ? 2 : 4); + + /*--- Initialize nVarGrad for deallocation ---*/ + + nVarGrad = nPrimVarGrad; + + nMarker = config->GetnMarker_All(); + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + /*--- Store the number of vertices on each marker for deallocation later ---*/ + + nVertex.resize(nMarker); + for (iMarker = 0; iMarker < nMarker; iMarker++) + nVertex[iMarker] = geometry->nVertex[iMarker]; + + /*--- Perform the non-dimensionalization for the flow equations using the + specified reference values. ---*/ + + // SetNondimensionalization(config, iMesh); + + /*--- Check if we are executing a verification case. If so, the + VerificationSolution object will be instantiated for a particular + option from the available library of verification solutions. Note + that this is done after SetNondim(), as problem-specific initial + parameters are needed by the solution constructors. ---*/ + + SetVerificationSolution(nDim, nVar, config); + + /*--- Allocate base class members. ---*/ + + Allocate(*config); + + /*--- MPI + OpenMP initialization. ---*/ + + HybridParallelInitialization(*config, *geometry); + + /*--- Jacobians and vector structures for implicit computations ---*/ + + if (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT) { + + if (rank == MASTER_NODE) + cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; + + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); + } + else { + if (rank == MASTER_NODE) + cout << "Explicit scheme. No Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; + } + + /*--- Read farfield conditions ---*/ + + Density_Inf = config->GetDensity_FreeStreamND(); + Pressure_Inf = config->GetPressure_FreeStreamND(); + Velocity_Inf = config->GetVelocity_FreeStreamND(); + Temperature_Inf = config->GetTemperature_FreeStreamND(); + if (config->GetKind_Species_Model() != SPECIES_MODEL::NONE) scalar_init = config->GetSpecies_Init(); + //GetFluidModel()->SetTDState_T(Temperature_Inf, scalar_init); + //Enthalpy_Inf = GetFluidModel()->GetEnthalpy(); + + /*--- Initialize the secondary values for direct derivative approximations ---*/ + + switch (config->GetDirectDiff()) { + case NO_DERIVATIVE: + /*--- Default ---*/ + break; + case D_DENSITY: + SU2_TYPE::SetDerivative(Density_Inf, 1.0); + break; + case D_PRESSURE: + SU2_TYPE::SetDerivative(Pressure_Inf, 1.0); + break; + case D_TEMPERATURE: + SU2_TYPE::SetDerivative(Temperature_Inf, 1.0); + break; + case D_MACH: case D_AOA: + case D_SIDESLIP: case D_REYNOLDS: + case D_TURB2LAM: case D_DESIGN: + /*--- Already done in postprocessing of config ---*/ + break; + default: + break; + } + + SetReferenceValues(*config); + + /*--- Initialize the solution to the far-field state everywhere. ---*/ + + if (navier_stokes) { + nodes = new CPBIncNSVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); + } else { + nodes = new CPBIncEulerVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); + } + SetBaseClassPointerToNodes(); + + if (iMesh == MESH_0) { + nodes->NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0; + } + + /*--- Initial comms. ---*/ + + CommunicateInitialState(geometry, config); + + /*--- Sizing edge mass flux array ---*/ + if (config->GetBounded_Scalar()) + EdgeMassFluxes.resize(geometry->GetnEdge()) = su2double(0.0); + + /*--- Add the solver name. ---*/ + SolverName = "INC.FLOW"; + + /*--- Finally, check that the static arrays will be large enough (keep this + * check at the bottom to make sure we consider the "final" values). ---*/ + if((nDim > MAXNDIM) || (nPrimVar > MAXNVAR)) + SU2_MPI::Error("Oops! The CPBIncEulerSolver static array sizes are not large enough.", CURRENT_FUNCTION); + +} diff --git a/SU2_CFD/src/solvers/CPBIncNSSolver.cpp b/SU2_CFD/src/solvers/CPBIncNSSolver.cpp new file mode 100644 index 00000000000..91d875acb15 --- /dev/null +++ b/SU2_CFD/src/solvers/CPBIncNSSolver.cpp @@ -0,0 +1,42 @@ +/*! + * \file CPBIncNSSolver.cpp + * \brief Main subroutines for solving pressure based Navier-Stokes incompressible flow. + * \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/CPBIncNSSolver.hpp" +#include "../../include/variables/CPBIncNSVariable.hpp" +#include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../include/solvers/CFVMFlowSolverBase.inl" + +/*--- Explicit instantiation of the parent class of CPBIncEulerSolver, + * to spread the compilation over two cpp files. ---*/ +template class CFVMFlowSolverBase; + + +CPBIncNSSolver::CPBIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : + CPBIncEulerSolver(geometry, config, iMesh, true) { + SU2_ZONE_SCOPED + +} diff --git a/SU2_CFD/src/solvers/CPoissonSolver.cpp b/SU2_CFD/src/solvers/CPoissonSolver.cpp new file mode 100644 index 00000000000..bc076bee896 --- /dev/null +++ b/SU2_CFD/src/solvers/CPoissonSolver.cpp @@ -0,0 +1,41 @@ +/*! + * \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) { + SU2_ZONE_SCOPED + +} + diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index a01590c6ab2..7899e5a2911 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -29,8 +29,11 @@ #include "../../include/solvers/CSolverFactory.hpp" #include "../../include/solvers/CEulerSolver.hpp" #include "../../include/solvers/CIncEulerSolver.hpp" +#include "../../include/solvers/CPBIncEulerSolver.hpp" #include "../../include/solvers/CNSSolver.hpp" #include "../../include/solvers/CIncNSSolver.hpp" +#include "../../include/solvers/CPBIncNSSolver.hpp" +#include "../../include/solvers/CPoissonSolver.hpp" #include "../../include/solvers/CNEMOEulerSolver.hpp" #include "../../include/solvers/CNEMONSSolver.hpp" #include "../../include/solvers/CTurbSASolver.hpp" @@ -69,6 +72,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 +85,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 +108,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 +150,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 +161,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 +174,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 +340,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() == ENUM_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; @@ -491,11 +505,31 @@ CSolver* CSolverFactory::CreateFlowSolver(SUB_SOLVER_TYPE kindFlowSolver, CSolve flowSolver = new CNSSolver(geometry, config, iMGLevel); break; case SUB_SOLVER_TYPE::INC_EULER: - flowSolver = new CIncEulerSolver(geometry, config, iMGLevel); + switch (config->GetKind_Incomp_System()) { + case ENUM_INCOMP_SYSTEM::DENSITY_BASED: + flowSolver = new CIncEulerSolver(geometry, config, iMGLevel); + break; + case ENUM_INCOMP_SYSTEM::PRESSURE_BASED: + flowSolver = new CPBIncEulerSolver(geometry, config, iMGLevel); + break; + default: + SU2_MPI::Error("Flow solver not found due to invalid incompressible method", CURRENT_FUNCTION); + break; + } flowSolver->Preprocessing(geometry, solver, config, iMGLevel, NO_RK_ITER, RUNTIME_FLOW_SYS, false); break; case SUB_SOLVER_TYPE::INC_NAVIER_STOKES: - flowSolver = new CIncNSSolver(geometry, config, iMGLevel); + switch (config->GetKind_Incomp_System()) { + case ENUM_INCOMP_SYSTEM::DENSITY_BASED: + flowSolver = new CIncNSSolver(geometry, config, iMGLevel); + break; + case ENUM_INCOMP_SYSTEM::PRESSURE_BASED: + flowSolver = new CPBIncNSSolver(geometry, config, iMGLevel); + break; + default: + SU2_MPI::Error("Flow solver not found due to invalid incompressible method", CURRENT_FUNCTION); + break; + } break; case SUB_SOLVER_TYPE::NEMO_EULER: flowSolver = new CNEMOEulerSolver(geometry, config, iMGLevel); diff --git a/SU2_CFD/src/variables/CPBIncEulerVariable.cpp b/SU2_CFD/src/variables/CPBIncEulerVariable.cpp new file mode 100644 index 00000000000..e38453ca623 --- /dev/null +++ b/SU2_CFD/src/variables/CPBIncEulerVariable.cpp @@ -0,0 +1,37 @@ +/*! + * \file CPBIncEulerVariable.cpp + * \brief Definition of the variable classes for pressure based incompressible flow. + * \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/CPBIncEulerVariable.hpp" +#include "../../include/fluid/CFluidModel.hpp" + +CPBIncEulerVariable::CPBIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy, + unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) + : CFlowVariable(npoint, ndim, nvar, ndim + 10, + ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config), + indices(ndim, 0) { + +} diff --git a/SU2_CFD/src/variables/CPBIncNSVariable.cpp b/SU2_CFD/src/variables/CPBIncNSVariable.cpp new file mode 100644 index 00000000000..26c66afbf60 --- /dev/null +++ b/SU2_CFD/src/variables/CPBIncNSVariable.cpp @@ -0,0 +1,38 @@ +/*! + * \file CPBIncNSVariable.cpp + * \brief Definition of the variable classes for incompressible flow. + * \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/CPBIncNSVariable.hpp" +#include "../../include/fluid/CFluidModel.hpp" + +CPBIncNSVariable::CPBIncNSVariable(su2double pressure, const su2double *velocity, su2double enthalpy, + unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) : + CPBIncEulerVariable(pressure, velocity, enthalpy, npoint, ndim, nvar, config), + Energy(config->GetEnergy_Equation()) { + + Vorticity.resize(nPoint, 3); + +} diff --git a/SU2_CFD/src/variables/CPoissonVariable.cpp b/SU2_CFD/src/variables/CPoissonVariable.cpp new file mode 100644 index 00000000000..e69de29bb2d