From df00fb7b7a34e08765796fd3342ac825e690694f Mon Sep 17 00:00:00 2001 From: bigboateng Date: Fri, 31 Jul 2026 10:06:19 +0100 Subject: [PATCH] Fix implicit target-time evaluation --- Common/include/CConfig.hpp | 30 +++++ Common/src/CConfig.cpp | 12 +- Common/src/grid_movement/CSurfaceMovement.cpp | 11 +- .../src/grid_movement/CVolumetricMovement.cpp | 17 +-- SU2_CFD/src/drivers/CDriver.cpp | 2 +- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 2 +- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 2 +- SU2_CFD/src/iteration/CFluidIteration.cpp | 12 +- SU2_CFD/src/output/COutput.cpp | 7 +- SU2_CFD/src/solvers/CMeshSolver.cpp | 20 ++- UnitTests/Common/CConfig_tests.cpp | 122 ++++++++++++++++++ UnitTests/meson.build | 1 + 12 files changed, 196 insertions(+), 42 deletions(-) create mode 100644 UnitTests/Common/CConfig_tests.cpp diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 1b2a45dc3ff..d5181f35bc1 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -3400,6 +3400,30 @@ class CConfig { */ su2double GetPhysicalTime(void) const { return PhysicalTime; } + /*! + * \brief Get the physical time represented by the completed state of the current time iteration. + * \return Nondimensional physical time at the end of the current time step. + */ + su2double GetPhysicalTimeAtEndOfTimeStep(void) const { + if (!ContinuousAdjoint && !DiscreteAdjoint && TimeMarching == TIME_MARCHING::TIME_STEPPING) + return PhysicalTime + Delta_UnstTimeND; + return PhysicalTime; + } + + /*! + * \brief Project a zero-based time iteration label onto the physical time represented by the solver state. + * \param[in] val_iter - Zero-based time iteration label used for output and restart files. + * \return Nondimensional physical time represented by the solver state. + */ + su2double GetPhysicalTime(unsigned long val_iter) const { + if (TimeMarching == TIME_MARCHING::STEADY) return 0.0; + if (!ContinuousAdjoint && !DiscreteAdjoint && + (TimeMarching == TIME_MARCHING::DT_STEPPING_1ST || + TimeMarching == TIME_MARCHING::DT_STEPPING_2ND)) + ++val_iter; + return static_cast(val_iter) * Delta_UnstTimeND; + } + /*! * \brief Get information about writing the performance summary at the end of a calculation. * \return TRUE means that the performance summary will be written at the end of a calculation. @@ -6146,6 +6170,12 @@ class CConfig { */ su2double GetCurrent_UnstTime(void) const { return Current_UnstTime; } + /*! + * \brief Get the current nondimensional physical time in an unsteady simulation. + * \return Value of the nondimensional physical time in an unsteady simulation. + */ + su2double GetCurrent_UnstTimeND(void) const { return Current_UnstTimeND; } + /*! * \brief Divide the rectbles and hexahedron. * \return TRUE if the elements must be divided; otherwise FALSE. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a9ca92ab145..5f4bba02275 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1138,6 +1138,9 @@ void CConfig::SetPointersNull() { InnerIter = 0; nIntCoeffs = 0; OuterIter = 0; + PhysicalTime = 0.0; + Current_UnstTime = 0.0; + Current_UnstTimeND = 0.0; AoA_Offset = 0; AoS_Offset = 0; @@ -8921,8 +8924,13 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, /*--- Set the simulation global time ---*/ - Current_UnstTime = static_cast(TimeIter)*Delta_UnstTime; - Current_UnstTimeND = static_cast(TimeIter)*Delta_UnstTimeND; + if (Time_Domain && !ContinuousAdjoint && !DiscreteAdjoint) { + Current_UnstTime = PhysicalTime * Time_Ref; + Current_UnstTimeND = PhysicalTime; + } else { + Current_UnstTime = static_cast(TimeIter)*Delta_UnstTime; + Current_UnstTimeND = static_cast(TimeIter)*Delta_UnstTimeND; + } /*--- Set the solver methods ---*/ diff --git a/Common/src/grid_movement/CSurfaceMovement.cpp b/Common/src/grid_movement/CSurfaceMovement.cpp index 6eb1b09ab40..eb4343bc120 100644 --- a/Common/src/grid_movement/CSurfaceMovement.cpp +++ b/Common/src/grid_movement/CSurfaceMovement.cpp @@ -3349,7 +3349,7 @@ void CSurfaceMovement::AeroelasticDeform(CGeometry* geometry, CConfig* config, u su2double Omega, dt, psi; dt = config->GetDelta_UnstTimeND(); Omega = config->GetRotation_Rate(2) / config->GetOmega_Ref(); - psi = Omega * (dt * TimeIter); + psi = Omega * config->GetPhysicalTime(); /*--- Correct for the airfoil starting position (This is hardcoded in here) ---*/ if (Monitoring_Tag == "Airfoil1") { @@ -3437,9 +3437,8 @@ void CSurfaceMovement::SetBoundary_Flutter3D(CGeometry* geometry, CConfig* confi } else { /*--- Forward time for the direct problem ---*/ - time_new = static_cast(iter) * deltaT; - time_old = time_new; - if (iter != 0) time_old = (static_cast(iter) - 1.0) * deltaT; + time_new = config->GetPhysicalTime(); + time_old = time_new > deltaT ? time_new - deltaT : 0.0; } /*--- Update the pitching angle at this time step. Flip sign for @@ -3614,8 +3613,8 @@ void CSurfaceMovement::SetExternal_Deformation(CGeometry* geometry, CConfig* con else dt = -1.0 * dt; } else { - /*--- No rotation at all for the first direct solution ---*/ - if (iter == 0) dt = 0; + /*--- Rotate the direct solution to the physical target time. ---*/ + dt = config->GetPhysicalTime(); } /*--- Compute delta change in the angle about the x, y, & z axes. ---*/ diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index ca6001b7f9f..63cc51ee5ba 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -507,7 +507,7 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry* geometry, CConfig* config, u dt = -1.0 * dt; } else { /*--- No rotation at all for the first direct solution ---*/ - if (iter == 0) dt = 0; + if (iter == 0 && config->GetPhysicalTime() == 0.0) dt = 0; } /*--- Center of rotation & angular velocity vector from config ---*/ @@ -674,13 +674,12 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry* geometry, CConfig* config, u if (iter != 0) time_old = (static_cast(directIter) + 1.0) * deltaT; } else { /*--- Forward time for the direct problem ---*/ - time_new = static_cast(iter) * deltaT; + time_new = harmonic_balance ? static_cast(iter) * deltaT : config->GetPhysicalTime(); if (harmonic_balance) { /*--- For harmonic balance, begin movement from the zero position ---*/ time_old = 0.0; } else { - time_old = time_new; - if (iter != 0) time_old = (static_cast(iter) - 1.0) * deltaT; + time_old = time_new > deltaT ? time_new - deltaT : 0.0; } } @@ -813,13 +812,12 @@ void CVolumetricMovement::Rigid_Plunging(CGeometry* geometry, CConfig* config, u if (iter != 0) time_old = (static_cast(directIter) + 1.0) * deltaT; } else { /*--- Forward time for the direct problem ---*/ - time_new = static_cast(iter) * deltaT; + time_new = harmonic_balance ? static_cast(iter) * deltaT : config->GetPhysicalTime(); if (harmonic_balance) { /*--- For harmonic balance, begin movement from the zero position ---*/ time_old = 0.0; } else { - time_old = time_new; - if (iter != 0) time_old = (static_cast(iter) - 1.0) * deltaT; + time_old = time_new > deltaT ? time_new - deltaT : 0.0; } } @@ -939,13 +937,12 @@ void CVolumetricMovement::Rigid_Translation(CGeometry* geometry, CConfig* config if (iter != 0) time_old = (static_cast(directIter) + 1.0) * deltaT; } else { /*--- Forward time for the direct problem ---*/ - time_new = static_cast(iter) * deltaT; + time_new = harmonic_balance ? static_cast(iter) * deltaT : config->GetPhysicalTime(); if (harmonic_balance) { /*--- For harmonic balance, begin movement from the zero position ---*/ time_old = 0.0; } else { - time_old = time_new; - if (iter != 0) time_old = (static_cast(iter) - 1.0) * deltaT; + time_old = time_new > deltaT ? time_new - deltaT : 0.0; } } diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 63d48590098..af637d7562e 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2956,7 +2956,7 @@ void CFluidDriver::Preprocess(unsigned long Iter) { for (iZone = 0; iZone < nZone; iZone++) { config_container[iZone]->SetInnerIter(Iter); if (config_container[iZone]->GetTime_Marching() != TIME_MARCHING::STEADY) - config_container[iZone]->SetPhysicalTime(static_cast(Iter)*config_container[iZone]->GetDelta_UnstTimeND()); + config_container[iZone]->SetPhysicalTime(config_container[iZone]->GetPhysicalTime(Iter)); else config_container[iZone]->SetPhysicalTime(0.0); } diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 00bdb364a44..c1ddf5fdff4 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -226,7 +226,7 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { general once the drivers are more stable. ---*/ if (driver_config->GetTime_Domain()) { - config_container[iZone]->SetPhysicalTime(static_cast(TimeIter)*config_container[iZone]->GetDelta_UnstTimeND()); + config_container[iZone]->SetPhysicalTime(config_container[iZone]->GetPhysicalTime(TimeIter)); } else { config_container[iZone]->SetPhysicalTime(0.0); diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index fa6e5ac3006..6e44c8c2597 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -117,7 +117,7 @@ void CSinglezoneDriver::Preprocess(unsigned long TimeIter) { general once the drivers are more stable. ---*/ if (config_container[ZONE_0]->GetTime_Marching() != TIME_MARCHING::STEADY) - config_container[ZONE_0]->SetPhysicalTime(static_cast(TimeIter)*config_container[ZONE_0]->GetDelta_UnstTimeND()); + config_container[ZONE_0]->SetPhysicalTime(config_container[ZONE_0]->GetPhysicalTime(TimeIter)); else config_container[ZONE_0]->SetPhysicalTime(0.0); diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index a79a9f1f8f4..e2d14022f24 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -498,11 +498,13 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C su2double x, y, x_gust, Gust[3] = {0.0}, NewGridVel[3] = {0.0}; const su2double* GridVel = nullptr; - su2double Physical_dt = config->GetDelta_UnstTime(); - unsigned long TimeIter = config->GetTimeIter(); - if (config->GetDiscrete_Adjoint()) TimeIter = config->GetUnst_AdjointIter() - TimeIter - 1; - - su2double Physical_t = TimeIter * Physical_dt; + su2double Physical_t; + if (config->GetDiscrete_Adjoint()) { + const auto TimeIter = config->GetUnst_AdjointIter() - config->GetTimeIter() - 1; + Physical_t = TimeIter * config->GetDelta_UnstTime(); + } else { + Physical_t = config->GetPhysicalTime() * config->GetTime_Ref(); + } su2double Uinf = solver[MESH_0][FLOW_SOL]->GetVelocity_Inf(0); // Assumption gust moves at infinity velocity diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index e210fbb59c4..001616c3429 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -2303,11 +2303,10 @@ void COutput::LoadCommonHistoryData(const CConfig *config) { SetHistoryOutputValue("TIME_STEP", config->GetDelta_UnstTimeND()*config->GetTime_Ref()); - /*--- Update the current time only if the time iteration has changed ---*/ + /*--- Report the absolute physical time represented by the completed state. + * Time iteration labels remain zero-based for output and restart files. ---*/ - if (SU2_TYPE::Int(GetHistoryFieldValue("TIME_ITER")) != static_cast(curTimeIter)) { - SetHistoryOutputValue("CUR_TIME", GetHistoryFieldValue("CUR_TIME") + GetHistoryFieldValue("TIME_STEP")); - } + SetHistoryOutputValue("CUR_TIME", config->GetPhysicalTimeAtEndOfTimeStep() * config->GetTime_Ref()); SetHistoryOutputValue("TIME_ITER", curTimeIter); SetHistoryOutputValue("INNER_ITER", curInnerIter); diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 35a633dc870..df27b164679 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -993,9 +993,8 @@ void CMeshSolver::Surface_Pitching(CGeometry *geometry, CConfig *config, unsigne /*--- Compute delta time based on physical time step ---*/ - time_new = iter*deltaT; - if (iter == 0) time_old = time_new; - else time_old = (iter-1)*deltaT; + time_new = config->GetPhysicalTime(); + time_old = time_new > deltaT ? time_new - deltaT : 0.0; /*--- Store displacement of each node on the pitching surface ---*/ /*--- Loop over markers and find the particular marker(s) (surface) to pitch ---*/ @@ -1126,9 +1125,8 @@ void CMeshSolver::Surface_Rotating(CGeometry *geometry, CConfig *config, unsigne /*--- Compute delta time based on physical time step ---*/ - time_new = iter*deltaT; - if (iter == 0) time_old = time_new; - else time_old = (iter-1)*deltaT; + time_new = config->GetPhysicalTime(); + time_old = time_new > deltaT ? time_new - deltaT : 0.0; /*--- Store displacement of each node on the rotating surface ---*/ /*--- Loop over markers and find the particular marker(s) (surface) to rotate ---*/ @@ -1307,9 +1305,8 @@ void CMeshSolver::Surface_Plunging(CGeometry *geometry, CConfig *config, unsigne /*--- Compute delta time based on physical time step ---*/ - time_new = iter*deltaT; - if (iter == 0) time_old = time_new; - else time_old = (iter-1)*deltaT; + time_new = config->GetPhysicalTime(); + time_old = time_new > deltaT ? time_new - deltaT : 0.0; /*--- Store displacement of each node on the plunging surface ---*/ /*--- Loop over markers and find the particular marker(s) (surface) to plunge ---*/ @@ -1434,9 +1431,8 @@ void CMeshSolver::Surface_Translating(CGeometry *geometry, CConfig *config, unsi /*--- Compute delta time based on physical time step ---*/ - time_new = iter*deltaT; - if (iter == 0) time_old = time_new; - else time_old = (iter-1)*deltaT; + time_new = config->GetPhysicalTime(); + time_old = time_new > deltaT ? time_new - deltaT : 0.0; /*--- Store displacement of each node on the translating surface ---*/ /*--- Loop over markers and find the particular marker(s) (surface) to translate ---*/ diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp new file mode 100644 index 00000000000..2dd1a2a9dac --- /dev/null +++ b/UnitTests/Common/CConfig_tests.cpp @@ -0,0 +1,122 @@ +/*! + * \file CConfig_tests.cpp + * \brief Unit tests for configuration-level time semantics. + * \author SU2 Contributors + * \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 + * + * 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 "catch.hpp" +#include +#include "../../Common/include/CConfig.hpp" +#include "../../Common/include/grid_movement/CVolumetricMovement.hpp" +#include "../UnitQuadTestCase.hpp" + +namespace { + +std::unique_ptr MakeTimeConfig(const char* time_marching, const char* math_problem = "DIRECT") { + std::stringstream options; + options << "SOLVER= EULER\n" + << "MATH_PROBLEM= " << math_problem << "\n" + << "TIME_DOMAIN= YES\n" + << "TIME_MARCHING= " << time_marching << "\n" + << "TIME_STEP= 0.5\n" + << "TIME_ITER= 10\n"; + return std::unique_ptr(new CConfig(options, SU2_COMPONENT::SU2_CFD, false)); +} + +} // namespace + +TEST_CASE("Physical time distinguishes a time-iteration label from its target state", "[CConfig][time]") { + constexpr su2double delta_time = 0.25; + + SECTION("dual-time labels denote the state at the end of the step") { + auto config = MakeTimeConfig("DUAL_TIME_STEPPING-1ST_ORDER"); + config->SetDelta_UnstTimeND(delta_time); + + CHECK(config->GetPhysicalTime(0) == Approx(delta_time)); + CHECK(config->GetPhysicalTime(7) == Approx(8.0 * delta_time)); + + config->SetTimeIter(0); + config->SetPhysicalTime(config->GetPhysicalTime(0)); + config->SetGlobalParam(config->GetKind_Solver(), RUNTIME_FLOW_SYS); + + CHECK(config->GetTimeIter() == 0); + CHECK(config->GetCurrent_UnstTimeND() == Approx(config->GetPhysicalTime())); + CHECK(config->GetPhysicalTimeAtEndOfTimeStep() == Approx(delta_time)); + } + + SECTION("second-order dual time has the same target-time convention") { + auto config = MakeTimeConfig("DUAL_TIME_STEPPING-2ND_ORDER"); + config->SetDelta_UnstTimeND(delta_time); + + CHECK(config->GetPhysicalTime(0) == Approx(delta_time)); + } + + SECTION("stage-based time stepping starts at the initial time") { + auto config = MakeTimeConfig("TIME_STEPPING"); + config->SetDelta_UnstTimeND(delta_time); + + CHECK(config->GetPhysicalTime(0) == Approx(0.0)); + CHECK(config->GetPhysicalTime(7) == Approx(7.0 * delta_time)); + + config->SetPhysicalTime(config->GetPhysicalTime(0)); + CHECK(config->GetPhysicalTimeAtEndOfTimeStep() == Approx(delta_time)); + } + + SECTION("continuous adjoint time retains its reverse-time interpretation") { + auto config = MakeTimeConfig("DUAL_TIME_STEPPING-1ST_ORDER", "CONTINUOUS_ADJOINT"); + config->SetDelta_UnstTimeND(delta_time); + + CHECK(config->GetPhysicalTime(0) == Approx(0.0)); + } +} + +TEST_CASE("Rigid motion evaluates the first dual-time target state", "[CConfig][time][grid movement]") { + UnitQuadTestCase test_case; + test_case.AddOption("TIME_DOMAIN= YES"); + test_case.AddOption("TIME_MARCHING= DUAL_TIME_STEPPING-1ST_ORDER"); + test_case.AddOption("TIME_STEP= 0.25"); + test_case.AddOption("TIME_ITER= 1"); + test_case.AddOption("GRID_MOVEMENT= RIGID_MOTION"); + test_case.AddOption("PLUNGING_OMEGA= 0.0 1.0 0.0"); + test_case.AddOption("PLUNGING_AMPL= 0.0 0.1 0.0"); + test_case.InitConfig(); + + constexpr su2double delta_time = 0.25; + test_case.config->SetDelta_UnstTimeND(delta_time); + test_case.config->SetLength_Ref(1.0); + test_case.config->SetOmega_Ref(1.0); + test_case.config->SetPhysicalTime(test_case.config->GetPhysicalTime(0)); + test_case.InitGeometry(); + + const auto initial_y = test_case.geometry->nodes->GetCoord(0, 1); + const auto amplitude = test_case.config->GetPlunging_Ampl(1) / test_case.config->GetLength_Ref(); + const auto omega = test_case.config->GetPlunging_Omega(1) / test_case.config->GetOmega_Ref(); + const auto expected_delta = -amplitude * sin(omega * delta_time); + + CVolumetricMovement movement(test_case.geometry.get()); + movement.Rigid_Plunging(test_case.geometry.get(), test_case.config.get(), ZONE_0, 0); + + CHECK(test_case.geometry->nodes->GetCoord(0, 1) - initial_y == Approx(expected_delta)); + CHECK(expected_delta != Approx(0.0)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index f8c22511b5f..4672c4d1408 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -7,6 +7,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'Common/geometry/dual_grid/CDualGrid_tests.cpp', 'Common/geometry/CGeometry_test.cpp', + 'Common/CConfig_tests.cpp', 'Common/toolboxes/CQuasiNewtonInvLeastSquares_tests.cpp', 'Common/toolboxes/C1DInterpolation_tests.cpp', 'Common/vectorization.cpp',