[JBCFD]

All,

We've been running SU2 against our commercial codes for a while now to build some confidence in the software before putting it forward as a suitable companion CFD package to our current commercial code.

As such, we've developed a couple of relevant test cases that we run when we update the code. We were a few versions behind on the code, so when we updated to v 5.0.0, we noticed a significant problem with one of our test cases and I was hoping I could get some thoughts on the genesis of the issue.

I've traced the issue back to a change between v. 4.1.0 (works) and v. 4.1.1 (problems).

The test case in question is an axisymmetric viscous supersonic converging/diverging nozzle with a mesh aft of the exit to capture the plume characteristics open to atmosphere. I've added the config file contents below.

In version 4.1.0, we set up the case to run with an adaptive CFL ranging from 1.0 to 5.75, with the decrement parameter set to 0 so that it can only ramp in the positive direction. This resulted in good convergence after approx. 10,000 to 12,000 iterations... which we deemed acceptable given the nature of the flow and the fact that we weren't using any convergence accelerators.

In version 4.1.1, the same setup results in non-physical points in the solution at about 600 iterations and divergence around 1000 iterations. We've played with the CFL parameters as well as a few others. It seems the newer versions of the code require us to run the problem at MUCH lower CFL's. Turning on the "VISCOUS_LIMITER" variables seems to "help" with higher CFL's... but still limits us to a CFL of around 2.0, which takes quite a big longer to converge (50,000 iterations approx.).

Any thoughts you might have as to what changes between the two versions may have resulted in this behavior would be greatly appreciated. We'd like to keep updating versions of the code and using them as the tool develops, but we need to understand the genesis of this issue before we can do that with confidence.

Thanks in advance,

JB


START INPUT FILE DATA:

% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
%
AXISYMMETRIC= YES
% Physical governing equations (EULER, NAVIER_STOKES,
% WAVE_EQUATION, HEAT_EQUATION, LINEAR_ELASTICITY,
% POISSON_EQUATION)
PHYSICAL_PROBLEM= NAVIER_STOKES
%
% Specify turbulence model (NONE, SA, SA_NEG, SST)
KIND_TURB_MODEL= SA
%
% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
MATH_PROBLEM= DIRECT
%
% Restart solution (NO, YES)
RESTART_SOL= NO
%
% Regime type (COMPRESSIBLE, INCOMPRESSIBLE)
REGIME_TYPE= COMPRESSIBLE
%
% System of measurements (SI, US)
% International system of units (SI): ( meters, kilograms, Kelvins,
% Newtons = kg m/s^2, Pascals = N/m^2,
% Density = kg/m^3, Speed = m/s,
% Equiv. Area = m^2 )
% United States customary units (US): ( inches, slug, Rankines, lbf = slug ft/s^2,
% psf = lbf/ft^2, Density = slug/ft^3,
% Speed = ft/s, Equiv. Area = ft^2 )
SYSTEM_MEASUREMENTS= SI
% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
%
% Mach number (non-dimensional, based on the free-stream values)
MACH_NUMBER= 0.01
%
% Angle of attack (degrees, only for compressible flows)
AoA= 0.0
%
% Free-stream option to choose between density and temperature (default) for
% initializing the solution (TEMPERATURE_FS, DENSITY_FS)
FREESTREAM_OPTION= TEMPERATURE_FS
%
% Free-stream pressure (101325.0 N/m^2, 2116.216 psf by default)
FREESTREAM_PRESSURE= 101325
%
% Free-stream temperature (288.15 K, 518.67 R by default)
FREESTREAM_TEMPERATURE= 300
%
% Reynolds number (non-dimensional, based on the free-stream values)
REYNOLDS_NUMBER= 1.87E6
%
% Reynolds length (1 m, 1 inch by default)
REYNOLDS_LENGTH= 1.0
% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
%
% Reference origin for moment computation (m or in)
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 (m or in)
REF_LENGTH_MOMENT= 1.0
%
% Reference area for force coefficients (0 implies automatic
% calculation) (m^2 or in^2)
REF_AREA= 1.0
%
% Flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE,
% FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE)
REF_DIMENSIONALIZATION= DIMENSIONAL
% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------%
%
% Different gas model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS)
FLUID_MODEL= IDEAL_GAS
%
% Ratio of specific heats (1.4 default and the value is hardcoded
% for the model STANDARD_AIR)
GAMMA_VALUE= 1.2311
%
% Specific gas constant (287.058 J/kg*K default and this value is hardcoded
% for the model STANDARD_AIR)
GAS_CONSTANT= 334
%
% Critical Temperature (131.00 K by default)
CRITICAL_TEMPERATURE= 131.00
%
% Critical Pressure (3588550.0 N/m^2 by default)
CRITICAL_PRESSURE= 3588550.0
%
% Critical Density (263.0 Kg/m3 by default)
CRITICAL_DENSITY= 263.0
%
% Acentri factor (0.035 (air))
ACENTRIC_FACTOR= 0.035
% --------------------------- VISCOSITY MODEL ---------------------------------%
%
% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY).
VISCOSITY_MODEL= SUTHERLAND
%
% Molecular Viscosity that would be constant (1.716E-5 by default)
MU_CONSTANT= 1.716E-5
%
% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI)
MU_REF= 8.4662E-6
%
% Sutherland Temperature Ref (273.15 K default value for AIR SI)
MU_T_REF= 2598
%
% Sutherland constant (110.4 default value for AIR SI)
SUTHERLAND_CONSTANT= 187
% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------%
%
% Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL).
CONDUCTIVITY_MODEL= CONSTANT_PRANDTL
%
% Molecular Thermal Conductivity that would be constant (0.0257 by default)
KT_CONSTANT= 0.0257

% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
%
% Navier-Stokes (no-slip), isothermal wall marker(s) (NONE = no marker)
% Format: ( marker name, constant wall temperature (K), ... )
MARKER_ISOTHERMAL= ( NozzleWall , 1000, Plenum, 300)
% Symmetry boundary marker(s) (NONE = no marker)
MARKER_SYM= ( Axis )
% Inlet boundary type (TOTAL_CONDITIONS, MASS_FLOW)
INLET_TYPE= TOTAL_CONDITIONS
%
% Inlet boundary marker(s) with the following formats (NONE = no marker)
% Total Conditions: (inlet marker, total temp, total pressure, flow_direction_x,
% flow_direction_y, flow_direction_z, ... ) where flow_direction is
% a unit vector.
% Mass Flow: (inlet marker, density, velocity magnitude, flow_direction_x,
% flow_direction_y, flow_direction_z, ... ) where flow_direction is
% a unit vector.
% Incompressible: (inlet marker, NULL, velocity magnitude, flow_direction_x,
% flow_direction_y, flow_direction_z, ... ) where flow_direction is
% a unit vector.
MARKER_INLET= ( Inlet, 2598, 8963184, 1, 0, 0, AirIn, 300, 101825, 0.707107, -0.707107, 0)
% Outlet boundary marker(s) (NONE = no marker)
% Format: ( outlet marker, back pressure (static), ... )
MARKER_OUTLET= (Outlet, 101325)

% ------------------------ SURFACES IDENTIFICATION ----------------------------%
%
% Marker(s) of the surface in the surface flow solution file
MARKER_PLOTTING = ( NozzleWall )
%
% Marker(s) of the surface where the non-dimensional coefficients are evaluated.
MARKER_MONITORING = ( NozzleWall )
%
% Marker(s) of the surface where obj. func. (design problem) will be evaluated
MARKER_DESIGNING = ( NozzleWall )
% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
%
% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
NUM_METHOD_GRAD= GREEN_GAUSS
%
% CFL number (stating value for the adaptive CFL number)
CFL_NUMBER= 1
%
% Adaptive CFL number (NO, YES)
CFL_ADAPT= YES
%
% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
% CFL max value )
CFL_ADAPT_PARAM= ( 0, 0.5, 0.9, 5.75 )
%
% Maximum Delta Time in local time stepping simulations
MAX_DELTA_TIME= 1E5
%
% Runge-Kutta alpha coefficients
RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
%
% Objective function in optimization problem (DRAG, LIFT, SIDEFORCE, MOMENT_X,
% MOMENT_Y, MOMENT_Z, EFFICIENCY,
% EQUIVALENT_AREA, NEARFIELD_PRESSURE,
% FORCE_X, FORCE_Y, FORCE_Z, THRUST,
% TORQUE, FREE_SURFACE, TOTAL_HEATFLUX,
% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE,
% INVERSE_DESIGN_HEATFLUX, AVG_TOTAL_PRESSURE,
% MASS_FLOW_RATE)

% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
%
% Reference element length for computing the slope and sharp edges
% limiters (0.1 m, 5.0 in by default)
REF_ELEM_LENGTH= 0.1
%
% Coefficient for the limiter
LIMITER_COEFF= 0.3
%
% Freeze the value of the limiter after a number of iterations
LIMITER_ITER= 999999
%
% Coefficient for the sharp edges limiter
SHARP_EDGES_COEFF= 3.0
%
% Reference coefficient (sensitivity) for detecting sharp edges.
REF_SHARP_EDGES= 3.0
%
% Remove sharp edges from the sensitivity evaluation (NO, YES)
SENS_REMOVE_SHARP= NO
% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
%
% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI,
% SMOOTHER_ILU0, SMOOTHER_LUSGS,
% SMOOTHER_LINELET)
LINEAR_SOLVER= FGMRES
%
% Preconditioner of the Krylov linear solver (ILU0, LU_SGS, LINELET, JACOBI)
LINEAR_SOLVER_PREC= LU_SGS
%
% Minimum error of the linear solver for implicit formulations
LINEAR_SOLVER_ERROR= 1E-4
%
% Max number of iterations of the linear solver for the implicit formulation
LINEAR_SOLVER_ITER= 5
% -------------------------- MULTIGRID PARAMETERS -----------------------------%
%
% Multi-grid levels (0 = no multi-grid)
MGLEVEL= 0
%
% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
MGCYCLE= W_CYCLE
%
% Multi-grid pre-smoothing level
MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
%
% Multi-grid post-smoothing level
MG_POST_SMOOTH= ( 0, 0, 0, 0 )
%
% Jacobi implicit smoothing of the correction
MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
%
% Damping factor for the residual restriction
MG_DAMP_RESTRICTION= 0.75
%
% Damping factor for the correction prolongation
MG_DAMP_PROLONGATION= 0.75
% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
%
% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
% TURKEL_PREC, MSW)
CONV_NUM_METHOD_FLOW= HLLC
%
% Spatial numerical order integration (1ST_ORDER, 2ND_ORDER, 2ND_ORDER_LIMITER)
SPATIAL_ORDER_FLOW= 2ND_ORDER_LIMITER
%
% Slope limiter (VENKATAKRISHNAN, BARTH_JESPERSEN)
SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
%
% Entropy fix coefficient (0.0 implies no entropy fixing, 1.0 implies scalar
% artificial dissipation)
ENTROPY_FIX_COEFF= 0.001
%
% 1st, 2nd and 4th order artificial dissipation coefficients
AD_COEFF_FLOW= ( 0.15, 0.5, 0.02 )
%
% Viscous limiter (NO, YES)
VISCOUS_LIMITER_FLOW= NO
%
% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
TIME_DISCRE_FLOW= EULER_IMPLICIT
%
% Relaxation coefficient
RELAXATION_FACTOR_FLOW= 1.0
% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
%
% Convective numerical method (SCALAR_UPWIND)
CONV_NUM_METHOD_TURB= SCALAR_UPWIND
%
% Spatial numerical order integration (1ST_ORDER, 2ND_ORDER, 2ND_ORDER_LIMITER)
SPATIAL_ORDER_TURB= 1ST_ORDER
%
% Slope limiter (VENKATAKRISHNAN)
SLOPE_LIMITER_TURB= VENKATAKRISHNAN
%
% Viscous limiter (NO, YES)
VISCOUS_LIMITER_TURB= NO
%
% Time discretization (EULER_IMPLICIT)
TIME_DISCRE_TURB= EULER_IMPLICIT
%
% Reduction factor of the CFL coefficient in the turbulence problem
CFL_REDUCTION_TURB= 1.0
%
% Relaxation coefficient
RELAXATION_FACTOR_TURB= 1.0
% --------------------- HEAT NUMERICAL METHOD DEFINITION ----------------------%
%
% Value of the thermal diffusivity
THERMAL_DIFFUSIVITY= 1.0


% --------------------------- CONVERGENCE PARAMETERS --------------------------%
%
% Number of total iterations
EXT_ITER= 25000
%
% Convergence criteria (CAUCHY, RESIDUAL)
%
CONV_CRITERIA= RESIDUAL
%
% Residual reduction (order of magnitude with respect to the initial value)
RESIDUAL_REDUCTION= 5
%
% Min value of the residual (log10 of the residual)
RESIDUAL_MINVAL= -8
%
% Start convergence criteria at iteration number
STARTCONV_ITER= 10
%
% Number of elements to apply the criteria
CAUCHY_ELEMS= 100
%
% Epsilon to control the series convergence
CAUCHY_EPS= 1E-10
%
% Direct function to apply the convergence criteria (LIFT, DRAG, NEARFIELD_PRESS)
CAUCHY_FUNC_FLOW= DRAG
%
% Adjoint function to apply the convergence criteria (SENS_GEOMETRY, SENS_MACH)
CAUCHY_FUNC_ADJFLOW= SENS_GEOMETRY
% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
%
% Mesh input file
MESH_FILENAME= ContouredSU2.su2
%
% Mesh input file format (SU2, CGNS)
MESH_FORMAT= SU2
%
% Mesh output file
MESH_OUT_FILENAME= mesh_out.su2
%
% Restart flow input file
SOLUTION_FLOW_FILENAME= solution_flow.dat
%
% Restart adjoint input file
SOLUTION_ADJ_FILENAME= solution_adj.dat
%
% Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW,
% FIELDVIEW, FIELDVIEW_BINARY)
OUTPUT_FORMAT= PARAVIEW
%
% Output file convergence history (w/o extension)
CONV_FILENAME= history
%
% Output file with the forces breakdown
BREAKDOWN_FILENAME= forces_breakdown.dat
%
% Output file restart flow
RESTART_FLOW_FILENAME= restart_flow.dat
%
% Output file restart adjoint
RESTART_ADJ_FILENAME= restart_adj.dat
%
% Output file flow (w/o extension) variables
VOLUME_FLOW_FILENAME= flow
%
% Output file adjoint (w/o extension) variables
VOLUME_ADJ_FILENAME= adjoint
%
% Output file surface flow coefficient (w/o extension)
SURFACE_FLOW_FILENAME= surface_flow
%
% Writing solution file frequency
WRT_SOL_FREQ= 1000
%
%
% Writing convergence history frequency
WRT_CON_FREQ= 1
%
% Output residual values in the solution files
WRT_RESIDUALS= NO
%
% Output limiters values in the solution files
WRT_LIMITERS= NO
%
% Output the sharp edges detector
WRT_SHARPEDGES= NO
%
% Minimize the required output memory
LOW_MEMORY_OUTPUT= NO
%
% Verbosity of console output: NONE removes minor MPI overhead (NONE, HIGH)
CONSOLE_OUTPUT_VERBOSITY= HIGH

[JBCFD]

Note: The HLLC method also seems to be involved. Switching methods to ROE tolerates a higher CFL number.

[Strat1994]

JB,

Unfortunately I cannot comment on the code changes that caused this issue. Rather I was wondering if you are willing or able to comment about your benchmarking efforts with SU2. I would like to undertake a similar effort and it would be nice to hear about your experience. Are there particular cases where SU2 excels or falls short of other codes in terms of accuracy or efficiency? How well does SU2 handle complex geometry and/or "bad' grids, etc? What parts of the code (solvers, settings, etc) have you tested most?

I realize that this is an open-ended question but any insight would be appreciated.

Thanks,
Andy

[JBCFD]

Andy,

Most of our problems involve high speed, high temperature flows, both internal and external, so these are the kinds of problems with which we've done the most benchmarking. Our approach thus far has been to develop a series of simple, small mesh test cases to run with each subsequent version of the code (to make sure no code changes have adversely affected the numerics/boundary conditions that we use) and then to run any real-problem mesh with both our commercial solver and SU2 and compare results.

To date, the biggest limitation is the lack of wall functions in SU2. Because of this, the wall grids must be much finer for boundary layer resolution and accurate calculation of associated values such as wall heat transfer. For bulk flow and determining things like plume structure, thrust etc SU2 seems to do just fine.

Also, the lack of multi-fluid and multi-phase capability can limit some of the work that we do here, but given that we're trying to solve real engineering problems at a fast pace, we rarely use this capability.

Finally, I'd say, don't expect the same quick convergence you get with a commercial code using all of their nifty convergence accelerators and all of that. To be honest, this isn't a big drain either, as there isn't that great of a difference between convergence times and SU2 is totally open source.

P.S. UPDATE TO THE ABOVE QUESTION I ASKED:
Seems both the Riemann solver and the exit boundary condition were giving me trouble. Making some modifications to both have seemingly fixed the issue. Moral of the story is, you may have to play around with numerics a bit to get the software to behave the way you want.

Cheers,
JB