From bffaaad34d982792cb710ba5baf1f690045d49ed Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 11:55:22 -0800 Subject: [PATCH 01/12] Add metric/adaptation config options and maps --- Common/include/CConfig.hpp | 115 ++++++++++++++++++++++++++++ Common/include/option_structure.hpp | 30 ++++++++ Common/src/CConfig.cpp | 105 +++++++++++++++++++++++++ 3 files changed, 250 insertions(+) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 794622bb1a1e..c2d7e854bd75 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1277,6 +1277,18 @@ class CConfig { /*--- Additional flamelet solver options ---*/ FluidFlamelet_ParsedOptions flamelet_ParsedOptions; /*!< \brief Additional flamelet solver options */ + /*--- Mesh adaptation options ---*/ + bool Compute_Metric; /*!< \brief Determines if error estimation is taking place */ + bool Normalize_Metric; /*!< \brief Determines if metric tensor normalization is taking place */ + unsigned short Kind_Hessian_Method; /*!< \brief Numerical method for computation of Hessians. */ + unsigned short nMetric_Sensor; /*!< \brief Number of sensors to use for adaptation. */ + METRIC_SENSOR* Metric_Sensor; /*!< \brief Sensors to use for adaptation. */ + unsigned short Metric_Norm; /*!< \brief Lp-norm for mesh adaptation */ + unsigned long Metric_Complexity; /*!< \brief Constraint mesh complexity */ + su2double Metric_Hmax, /*!< \brief Maximum cell size */ + Metric_Hmin, /*!< \brief Minimum cell size */ + Metric_ARmax; /*!< \brief Maximum cell aspect ratio */ + /*! * \brief Set the default values of config options not set in the config file using another config object. * \param config - Config object to use the default values from. @@ -10142,4 +10154,107 @@ class CConfig { */ const FluidFlamelet_ParsedOptions& GetFlameletParsedOptions() const { return flamelet_ParsedOptions; } + /*! + * \brief Check if a metric tensor field is being computed + * \return TRUE<\code> if a metric tensor field is being computed + */ + bool GetCompute_Metric(void) const { return Compute_Metric; } + + /*! + * \brief Check if metric tensor normalization is being carried out + * \return TRUE<\code> if metric normalization is taking place + */ + bool GetNormalize_Metric(void) const { return Normalize_Metric; } + + /*! + * \brief Check if goal-oriented error estimation is being carried out + * \return TRUE<\code> if goal-oriented error estimation is taking place + */ + bool GetGoal_Oriented_Metric(void) const { return (nMetric_Sensor > 0) && (Metric_Sensor[0] == METRIC_SENSOR::GOAL); } + + /*! + * \brief Get the kind of method for computation of Hessians used for anisotropy. + * \return Numerical method for computation of Hessians used for anisotropy. + */ + unsigned short GetKind_Hessian_Method(void) const { return Kind_Hessian_Method; } + + /*! + * \brief Get adaptation sensor + */ + METRIC_SENSOR GetMetric_Sensor(unsigned short iSens) const { return Metric_Sensor[iSens]; } + + /*! + * \brief Get corresponding string from metric sensor type + */ + string GetMetric_SensorString(unsigned short iSens) const { + string sensor_name; + switch (Metric_Sensor[iSens]) { + case METRIC_SENSOR::DENSITY: + sensor_name = "Density"; + break; + case METRIC_SENSOR::MACH: + sensor_name = "Mach"; + break; + case METRIC_SENSOR::PRESSURE: + sensor_name = "Pressure"; + break; + case METRIC_SENSOR::TOTAL_PRESSURE: + sensor_name = "Total Pressure"; + break; + case METRIC_SENSOR::TEMPERATURE: + sensor_name = "Temperature"; + break; + case METRIC_SENSOR::TEMPERATURE_VE: + sensor_name = "Temperature_ve"; + break; + case METRIC_SENSOR::ENERGY: + sensor_name = "Energy"; + break; + case METRIC_SENSOR::ENERGY_VE: + sensor_name = "Energy_ve"; + break; + case METRIC_SENSOR::GOAL: + sensor_name = "Goal"; + break; + default: + SU2_MPI::Error("Unsupported metric sensor.", CURRENT_FUNCTION); + } + + return sensor_name; + } + + /*! + * \brief Get number of adaptation sensors + */ + unsigned short GetnMetric_Sensor(void) const { return nMetric_Sensor; } + + /*! + * \brief Get adaptation norm value (Lp) + */ + unsigned short GetMetric_Norm(void) const { return Metric_Norm; } + + /*! + * \brief Get maximum cell size + * \return Maximum cell size + */ + su2double GetMetric_Hmax(void) const { return Metric_Hmax; } + + /*! + * \brief Get minimum cell size + * \return Minimum cell size + */ + su2double GetMetric_Hmin(void) const { return Metric_Hmin; } + + /*! + * \brief Get maximum cell aspect ratio + * \return Maximum cell aspect ratio + */ + su2double GetMetric_ARmax(void) const { return Metric_ARmax; } + + /*! + * \brief Get constraint complexity + * \return Mesh complexity + */ + unsigned long GetMetric_Complexity(void) const { return Metric_Complexity; } + }; diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b1387..6b7320b62bed 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2621,6 +2621,8 @@ 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_GRAD_ADAPT , /*!< \brief Gradient vectors for anisotropic sizing metric (periodic only). */ + PERIODIC_HESSIAN , /*!< \brief Hessian tensors for anisotropic sizing metric (periodic only). */ }; /*! @@ -2652,6 +2654,8 @@ 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. */ + GRADIENT_ADAPT , /*!< \brief Gradient vectors for anisotropic sizing metric tensor. */ + HESSIAN , /*!< \brief Hessian tensors for anisotropic sizing metric tensor. */ }; /*! @@ -2797,6 +2801,32 @@ static const MapType Sobolev_Modus_Map = { MakePair("ONLY_GRADIENT", ENUM_SOBOLEV_MODUS::ONLY_GRAD) }; +/*! + * \brief Types of sensors for anisotropic metric + */ +enum class METRIC_SENSOR { + DENSITY, /*!< \brief Density feature-based metric. */ + MACH, /*!< \brief Mach feature-based metric. */ + PRESSURE, /*!< \brief Pressure feature-based metric. */ + TOTAL_PRESSURE, /*!< \brief Total pressure feature-based metric. */ + TEMPERATURE, /*!< \brief Temperature feature-based metric. */ + TEMPERATURE_VE, /*!< \brief Vibrational/electronic temperature feature-based metric. */ + ENERGY, /*!< \brief Energy feature-based metric. */ + ENERGY_VE, /*!< \brief Vibrational/electronic energy feature-based metric. */ + GOAL, /*!< \brief Goal-oriented metric. */ +}; +static const MapType Metric_Sensor_Map = { + MakePair("DENSITY", METRIC_SENSOR::DENSITY) + MakePair("MACH", METRIC_SENSOR::MACH) + MakePair("PRESSURE", METRIC_SENSOR::PRESSURE) + MakePair("TOTAL_PRESSURE", METRIC_SENSOR::TOTAL_PRESSURE) + MakePair("TEMPERATURE", METRIC_SENSOR::TEMPERATURE) + MakePair("TEMPERATURE_VE", METRIC_SENSOR::TEMPERATURE_VE) + MakePair("ENERGY", METRIC_SENSOR::ENERGY) + MakePair("ENERGY_VE", METRIC_SENSOR::ENERGY_VE) + MakePair("GOAL", METRIC_SENSOR::GOAL) +}; + /*! * \brief Type of mesh deformation */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 0a9cca0f63f7..5232edc31896 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3034,6 +3034,55 @@ void CConfig::SetConfig_Options() { /*!\brief ROM_SAVE_FREQ \n DESCRIPTION: How often to save snapshots for unsteady problems.*/ addUnsignedShortOption("ROM_SAVE_FREQ", rom_save_freq, 1); + /*--- options that are used for mesh adaptation ---*/ + /*!\par CONFIG_CATEGORY:Adaptation Options \ingroup Config*/ + + /*!\brief COMPUTE_METRIC \n DESCRIPTION: Compute a metric tensor field */ + addBoolOption("COMPUTE_METRIC", Compute_Metric, false); + /*!\brief NORMALIZE_METRIC \n DESCRIPTION: Normalize the metric tensor */ + addBoolOption("NORMALIZE_METRIC", Normalize_Metric, true); + /*!\brief NUM_METHOD_HESS \n DESCRIPTION: Numerical method for Hessian computation \n OPTIONS: See \link Gradient_Map \endlink. \n DEFAULT: GREEN_GAUSS. \ingroup Config*/ + addEnumOption("NUM_METHOD_HESS", Kind_Hessian_Method, Gradient_Map, GREEN_GAUSS); + + /*!\brief METRIC_SENSOR \n DESCRIPTION: Sensors for mesh adaptation */ + addEnumListOption("METRIC_SENSOR", nMetric_Sensor, Metric_Sensor, Metric_Sensor_Map); + /*!\brief METRIC_NORM \n DESCRIPTION: Lp-norm for mesh adaptation */ + addUnsignedShortOption("METRIC_NORM", Metric_Norm, 2); + /*!\brief METRIC_COMPLEXITY \n DESCRIPTION: Constraint mesh complexity */ + addUnsignedLongOption("METRIC_COMPLEXITY", Metric_Complexity, 10000); + + /*!\brief METRIC_HMAX \n DESCRIPTION: Constraint maximum cell size */ + addDoubleOption("METRIC_HMAX", Metric_Hmax, 10.0); + /*!\brief METRIC_HMIN \n DESCRIPTION: Constraint minimum cell size */ + addDoubleOption("METRIC_HMIN", Metric_Hmin, 1.0E-8); + /*!\brief METRIC_ARMAX \n DESCRIPTION: Constraint maximum cell aspect ratio */ + addDoubleOption("METRIC_ARMAX", Metric_ARmax, 1.0E6); + + /*!\brief ADAP_ITER \n DESCRIPTION: Mesh adaptation inner iterations per complexity */ + addPythonOption("ADAP_ITER"); + /*!\brief ADAP_COMPLEXITIES \n DESCRIPTION: List of constraint (target) mesh complexities for mesh convergence study */ + addPythonOption("ADAP_COMPLEXITIES"); + /*!\brief ADAP_FLOW_ITERS \n DESCRIPTION: Primal solver iterations at each target complexity */ + addPythonOption("ADAP_FLOW_ITERS"); + /*!\brief ADAP_ADJ_ITERS \n DESCRIPTION: Adjoint solver iterations at each target complexity */ + addPythonOption("ADAP_ADJ_ITERS"); + /*!\brief ADAP_FLOW_CFLS \n DESCRIPTION: Primal solver CFL number at each target complexity */ + addPythonOption("ADAP_FLOW_CFLS"); + /*!\brief ADAP_ADJ_CFLS \n DESCRIPTION: Adjoint solver CFL number at each target complexity */ + addPythonOption("ADAP_ADJ_CFLS"); + /*!\brief ADAP_RESIDUAL_REDUCTIONS \n DESCRIPTION: Residual reduction at each target complexity */ + addPythonOption("ADAP_RESIDUAL_REDUCTIONS"); + /*!\brief ADAP_ITER \n DESCRIPTION: Maximum cell size at each target complexity */ + addPythonOption("ADAP_HMAXS"); + /*!\brief ADAP_ITER \n DESCRIPTION: Minimum cell size at each target complexity */ + addPythonOption("ADAP_HMINS"); + /*!\brief ADAP_HGRAD \n DESCRIPTION: Size gradation smoothing parameter */ + addPythonOption("ADAP_HGRAD"); + /*!\brief ADAP_ITER \n DESCRIPTION: Hausdorff distance parameter for surface remeshing */ + addPythonOption("ADAP_HAUSD"); + /*!\brief ADAP_ITER \n DESCRIPTION: A mesh adaptation option */ + addPythonOption("ADAP_ANGLE"); + /* END_CONFIG_OPTIONS */ } @@ -5752,6 +5801,38 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error("BOUNDED_SCALAR discretization can only be used for incompressible problems.", CURRENT_FUNCTION); } + /*--- Checks for mesh adaptation ---*/ + if (Compute_Metric) { + /*--- Check that config is valid for requested sensor ---*/ + for (auto iSensor = 0; iSensor < nMetric_Sensor; iSensor++) { + /*--- If using GOAL, it must be the only sensor and the discrete adjoint must be used ---*/ + if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL) { + if (nMetric_Sensor != 1) + SU2_MPI::Error("Adaptation sensor GOAL cannot be used with other sensors.", CURRENT_FUNCTION); + if (!DiscreteAdjoint) + SU2_MPI::Error("Adaptation sensor GOAL can only be computed for MATH_PROBLEM = DISCRETE_ADJOINT.", CURRENT_FUNCTION); + } + + if (Kind_Solver == MAIN_SOLVER::NEMO_EULER || Kind_Solver == MAIN_SOLVER::NEMO_NAVIER_STOKES) { + if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY) + SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for NEMO problems."), CURRENT_FUNCTION); + } + if (Kind_Solver != MAIN_SOLVER::NEMO_EULER && Kind_Solver != MAIN_SOLVER::NEMO_NAVIER_STOKES) { + if (Metric_Sensor[iSensor] == METRIC_SENSOR::TEMPERATURE_VE || Metric_Sensor[iSensor] == METRIC_SENSOR::ENERGY_VE) + SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for non-NEMO problems."), CURRENT_FUNCTION); + } + if (Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) { + if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::ENERGY) + SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for INC problems."), CURRENT_FUNCTION); + } + } + + /*--- Only GG Hessians for now ---*/ + if (Kind_Hessian_Method != GREEN_GAUSS) { + SU2_MPI::Error("NUM_METHOD_HESS must be GREEN_GAUSS.", CURRENT_FUNCTION); + } + } + } void CConfig::SetMarkers(SU2_COMPONENT val_software) { @@ -7925,6 +8006,30 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "Actuator disk BEM method propeller data read from file: " << GetBEM_prop_filename() << endl; } } + + if (val_software == SU2_COMPONENT::SU2_CFD || val_software == SU2_COMPONENT::SU2_SOL) { + if (Compute_Metric) { + cout << endl <<"---------------- Mesh Adaptation Information ( Zone " << iZone << " ) -----------------" << endl; + cout << "Adaptation sensor(s): "; + for (auto iSensor = 0; iSensor < nMetric_Sensor; iSensor++) { + cout << GetMetric_SensorString(iSensor); + if (iSensor < nMetric_Sensor - 1 ) cout << ", "; + } + cout << endl; + switch (Kind_Hessian_Method) { + case GREEN_GAUSS: cout << "Hessian for adaptive metric: Green-Gauss." << endl; break; + } + if (Normalize_Metric) { + cout << "Target complexity: " << Metric_Complexity << endl; + cout << "Lp norm: " << Metric_Norm << endl; + cout << "Min. edge length: " << Metric_Hmin << endl; + cout << "Max. edge length: " << Metric_Hmax << endl; + } + else { + cout << "Output unnormalized metric field." << endl; + } + } + } } bool CConfig::TokenizeString(string & str, string & option_name, From 3f2a7ddfd182e9cc314a7e32c7687cd4339f2261 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 11:59:11 -0800 Subject: [PATCH 02/12] Add adaptation primitive/gradient/Hessian and metric tensor to CVariable with initialization --- SU2_CFD/include/variables/CVariable.hpp | 118 ++++++++++++++++++++++++ SU2_CFD/src/variables/CFlowVariable.cpp | 7 ++ SU2_CFD/src/variables/CVariable.cpp | 8 ++ 3 files changed, 133 insertions(+) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index dc905d79fcb4..5cd23008e271 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -103,6 +103,11 @@ class CVariable { VectorType SolutionExtra_BGS_k; /*!< \brief Intermediate storage, enables cross term extraction as that is also pushed to Solution. */ + MatrixType Primitive_Adapt; /*!< \brief Variables for which we need gradients for anisotropy in mesh adaptation. */ + CVectorOfMatrix Gradient_Adapt; /*!< \brief Gradient of sensor used for anisotropy in mesh adaptation. */ + CVectorOfMatrix Hessian; /*!< \brief Hessian of sensor used for anisotropy in mesh adaptation. */ + su2matrix Metric; /*!< \brief Metric tensor used for anisotropy in mesh adaptation. */ + protected: unsigned long nPoint = 0; /*!< \brief Number of points in the domain. */ unsigned long nDim = 0; /*!< \brief Number of dimension of the problem. */ @@ -112,6 +117,8 @@ class CVariable { unsigned long nSecondaryVar = 0; /*!< \brief Number of secondary variables. */ unsigned long nSecondaryVarGrad = 0; /*!< \brief Number of secondaries for which a gradient is computed. */ unsigned long nAuxVar = 0; /*!< \brief Number of auxiliary variables. */ + unsigned long nSymMat = 0; /*!< \brief Number of symmetric matrix componenents for Hessian and metric tensor. */ + /*--- Only allow default construction by derived classes. ---*/ CVariable() = default; @@ -2371,4 +2378,115 @@ class CVariable { inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; } inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; } + + /*! + * \brief Set the gradient of the solution. + * \param[in] iPoint - Point index. + * \param[in] gradient - Gradient of the solution. + */ + inline void SetPrimitive_Adapt(unsigned long iPoint, unsigned long iVar, su2double primitive) { Primitive_Adapt(iPoint,iVar) = primitive; } + + /*! + * \brief Get the gradient of the entire solution. + * \return Reference to gradient. + */ + inline const MatrixType& GetPrimitive_Adapt(void) const { return Primitive_Adapt; } + + /*! + * \brief Get the value of the solution gradient. + * \param[in] iPoint - Point index. + * \return Value of the gradient solution. + */ + inline su2double *GetPrimitive_Adapt(unsigned long iPoint) { return Primitive_Adapt[iPoint]; } + + /*! + * \brief Get the value of the solution gradient. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \return Value of the solution gradient. + */ + inline su2double GetPrimitive_Adapt(unsigned long iPoint, unsigned long iVar) const { return Primitive_Adapt(iPoint,iVar); } + + /*! + * \brief Set the gradient of the solution. + * \param[in] iPoint - Point index. + * \param[in] gradient - Gradient of the solution. + */ + inline void SetGradient_Adapt(unsigned long iPoint, su2double** gradient) { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + for (unsigned long iDim = 0; iDim < nDim; iDim++) + Gradient_Adapt(iPoint,iVar,iDim) = gradient[iVar][iDim]; + } + + /*! + * \brief Get the gradient of the entire solution. + * \return Reference to gradient. + */ + inline CVectorOfMatrix& GetGradient_Adapt(void) { return Gradient_Adapt; } + + /*! + * \brief Get the value of the solution gradient. + * \param[in] iPoint - Point index. + * \return Value of the gradient solution. + */ + inline CMatrixView GetGradient_Adapt(unsigned long iPoint) { return Gradient_Adapt[iPoint]; } + + /*! + * \brief Get the value of the solution gradient. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \param[in] iDim - Dimension index. + * \return Value of the solution gradient. + */ + inline su2double GetGradient_Adapt(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Gradient_Adapt(iPoint,iVar,iDim); } + + /*! + * \brief Set the Hessian of the solution. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \param[in] iMat - Hessian tensor index. + * \param[in] hess - Hessian of the solution. + */ + inline void SetHessian(unsigned long iPoint, unsigned long iVar, unsigned long iMat, su2double hess) { Hessian(iPoint,iVar,iMat) = hess; } + + /*! + * \brief Get the Hessian tensor field. + * \return Reference to the Hessian field. + */ + inline CVectorOfMatrix& GetHessian(void) { return Hessian; } + + /*! + * \brief Get the value of the Hessian. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \param[in] iMat - Hessian tensor index. + * \return Value of the Hessian. + */ + inline su2double GetHessian(unsigned long iPoint, unsigned long iVar, unsigned long iMat) const { return Hessian(iPoint,iVar,iMat); } + + /*! + * \brief Set the value of the metric. + * \param[in] iMat - Metric tensor index. + * \param[in] metric - Metric value. + */ + inline void SetMetric(unsigned long iPoint, unsigned short iMat, double metric) { Metric(iPoint,iMat) = metric; } + + /*! + * \brief Get the metric of the entire solution. + * \return Reference to the metric tensor field. + */ + inline su2matrix& GetMetric(void) { return Metric; } + + /*! + * \brief Add the value of the metric. + * \param[in] iMat - Metric tensor index. + * \param[in] metric - Metric value. + */ + inline void AddMetric(unsigned long iPoint, unsigned short iMat, double metric) { Metric(iPoint,iMat) += metric; } + + /*! + * \brief Get the value of the metric. + * \param[in] iMat - Metric tensor index. + */ + inline double GetMetric(unsigned long iPoint, unsigned short iMat) const { return Metric(iPoint,iMat); } }; diff --git a/SU2_CFD/src/variables/CFlowVariable.cpp b/SU2_CFD/src/variables/CFlowVariable.cpp index 6d07949d79f2..d9350462ae19 100644 --- a/SU2_CFD/src/variables/CFlowVariable.cpp +++ b/SU2_CFD/src/variables/CFlowVariable.cpp @@ -97,6 +97,13 @@ CFlowVariable::CFlowVariable(unsigned long npoint, unsigned long ndim, unsigned if (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE) { HB_Source.resize(nPoint, nVar) = su2double(0.0); } + + if (config->GetCompute_Metric()) { + unsigned short nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + unsigned short nSymMat = 3 * (nDim - 1); + Primitive_Adapt.resize(nPoint, nHess) = su2double(0.0); + Metric.resize(nPoint, nSymMat) = 0.0; + } } void CFlowVariable::SetSolution_New() { diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 111a2a17775d..8e7234b258ca 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -51,6 +51,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva nPoint = npoint; nDim = ndim; nVar = nvar; + nSymMat = 3 * (nDim - 1); /*--- Allocate fields common to all problems. Do not allocate fields that are specific to one solver, i.e. not common, in this class. ---*/ @@ -79,6 +80,13 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if (config->GetMultizone_Problem()) Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0); + + /*--- Gradient and Hessian for anisotropic metric ---*/ + if (config->GetCompute_Metric()) { + unsigned short nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + Gradient_Adapt.resize(nPoint,nHess,nDim,0.0); + Hessian.resize(nPoint,nHess,nSymMat,0.0); + } } void CVariable::Set_OldSolution() { From 10e16068d5776f888df6700516e065e9c3b8bcfc Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 12:03:02 -0800 Subject: [PATCH 03/12] Add solver functionality to calculate metric tensor field - GG gradient/Hessian calculation on adaptation-specific variable containers - Metric integration and normalization for steady and unsteady flows - Storing of proper primitive variables for Euler solver --- .../gradients/computeGradientsGreenGauss.hpp | 137 +++++++ SU2_CFD/include/gradients/computeMetrics.hpp | 377 ++++++++++++++++++ SU2_CFD/include/solvers/CEulerSolver.hpp | 43 ++ SU2_CFD/include/solvers/CSolver.hpp | 63 ++- SU2_CFD/src/solvers/CBaselineSolver.cpp | 1 + SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 1 + SU2_CFD/src/solvers/CEulerSolver.cpp | 1 + SU2_CFD/src/solvers/CIncEulerSolver.cpp | 1 + SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 1 + SU2_CFD/src/solvers/CSolver.cpp | 133 +++++- SU2_CFD/src/solvers/CTurbSASolver.cpp | 3 +- 11 files changed, 744 insertions(+), 17 deletions(-) create mode 100644 SU2_CFD/include/gradients/computeMetrics.hpp diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp index 98a170cf3810..a2f984e1560a 100644 --- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp +++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp @@ -179,6 +179,123 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER solver->InitiateComms(&geometry, &config, kindMpiComm); solver->CompleteComms(&geometry, &config, kindMpiComm); } + +template +void computeHessiansGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PERIODIC_QUANTITIES kindPeriodicComm, + CGeometry& geometry, const CConfig& config, const GradientType& gradient, + const size_t varBegin, const size_t varEnd, const int idxVel, GradientType& hessian) { + const size_t nPointDomain = geometry.GetnPointDomain(); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + const auto chunkSize = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + + const size_t nSymMat = 3 * (nDim - 1); + su2double diagScale; + + /*--- For each (non-halo) volume integrate over its faces (edges). ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) { + auto nodes = geometry.nodes; + + /*--- Clear the Hessian. --*/ + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iMat = 0; iMat < nSymMat; ++iMat) hessian(iPoint, iVar, iMat) = 0.0; + + /*--- Handle averaging and division by volume in one constant. ---*/ + + su2double halfOnVol = 0.5 / (nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint)); + + /*--- Add a contribution due to each neighbor. ---*/ + + for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) { + size_t iEdge = nodes->GetEdge(iPoint,iNeigh); + size_t jPoint = nodes->GetPoint(iPoint,iNeigh); + + /*--- Determine if edge points inwards or outwards of iPoint. + * If inwards we need to flip the area vector. ---*/ + + su2double dir = (iPoint < jPoint)? 1.0 : -1.0; + su2double weight = dir * halfOnVol; + + const auto area = geometry.edges->GetNormal(iEdge); + + for (size_t jDim = 0; jDim < nDim; ++jDim) { + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) { + su2double flux = weight * (gradient(iPoint, iVar, jDim) + gradient(jPoint, iVar, jDim)); + for (size_t iDim = 0; iDim < nDim; ++iDim) { + size_t ind = (iDim <= jDim) ? iDim*nDim - ((iDim - 1)*iDim)/2 + jDim - iDim + : jDim*nDim - ((jDim - 1)*jDim)/2 + iDim - jDim; + diagScale = (iDim == jDim)? 1.0 : 0.5; + hessian(iPoint, iVar, ind) += diagScale * flux * area[iDim]; + } // iDims + } // variables + } // jDims + } // neighbors + } // points + END_SU2_OMP_FOR + + /*--- Add edges of markers that contribute to the Hessians ---*/ + for (size_t iMarker = 0; iMarker < geometry.GetnMarker(); ++iMarker) { + if ((config.GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + (config.GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) && + (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + /*--- Work is shared in inner loop as two markers + * may try to update the same point. ---*/ + + SU2_OMP_FOR_STAT(32) + for (size_t iVertex = 0; iVertex < geometry.GetnVertex(iMarker); ++iVertex) { + size_t iPoint = geometry.vertex[iMarker][iVertex]->GetNode(); + auto nodes = geometry.nodes; + + /*--- Halo points do not need to be considered. ---*/ + + if (!nodes->GetDomain(iPoint)) continue; + + su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint); + const auto area = geometry.vertex[iMarker][iVertex]->GetNormal(); + + for (size_t jDim = 0; jDim < nDim; ++jDim) { + for (size_t iVar = varBegin; iVar < varEnd; iVar++) { + su2double flux = gradient(iPoint, iVar, jDim) / volume; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + size_t ind = (iDim <= jDim) ? iDim * nDim - ((iDim - 1) * iDim)/2 + jDim - iDim + : jDim * nDim - ((jDim - 1) * jDim)/2 + iDim - jDim; + diagScale = (iDim == jDim)? 1.0 : 0.5; + hessian(iPoint, iVar, ind) -= diagScale * flux * area[iDim]; + } // iDims + } // variables + } // jDims + } // vertices + END_SU2_OMP_FOR + } //found right marker + } // iMarkers + + /*--- Compute the corrections for symmetry planes and Euler walls. ---*/ + /*--- TODO? correctHessiansSymmetry() ---*/ + + /*--- If no solver was provided we do not communicate ---*/ + + if (solver == nullptr) return; + + /*--- Account for periodic contributions. ---*/ + + for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) + { + solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + } + + /*--- Obtain the gradients at halo points from the MPI ranks that own them. ---*/ + + solver->InitiateComms(&geometry, &config, kindMpiComm); + solver->CompleteComms(&geometry, &config, kindMpiComm); + +} } // namespace detail @@ -205,3 +322,23 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER break; } } + +template +void computeHessiansGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PERIODIC_QUANTITIES kindPeriodicComm, + CGeometry& geometry, const CConfig& config, const GradientType& gradient, + const size_t varBegin, const size_t varEnd, const int idxVel, GradientType& hessian) { + switch (geometry.GetnDim()) { + case 2: + detail::computeHessiansGreenGauss<2>(solver, kindMpiComm, kindPeriodicComm, geometry, config, gradient, + varBegin, varEnd, idxVel, hessian); + break; + case 3: + detail::computeHessiansGreenGauss<3>(solver, kindMpiComm, kindPeriodicComm, geometry, config, gradient, + varBegin, varEnd, idxVel, hessian); + break; + default: + SU2_MPI::Error("Too many dimensions to compute Hessians.", CURRENT_FUNCTION); + break; + } +} + diff --git a/SU2_CFD/include/gradients/computeMetrics.hpp b/SU2_CFD/include/gradients/computeMetrics.hpp new file mode 100644 index 000000000000..cb437fbd29ea --- /dev/null +++ b/SU2_CFD/include/gradients/computeMetrics.hpp @@ -0,0 +1,377 @@ +/*! + * \file computeMetrics.hpp + * \brief Generic implementation of the metric tensor computation. + * \note This allows the same implementation to be used for goal-oriented + * or feature-based mesh adaptation. + * \author B. Munguía + * \version 8.3.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2025, 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 +#include +#include "../../../Common/include/parallelization/omp_structure.hpp" +#include "../../../Common/include/linear_algebra/blas_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +namespace tensor { + struct metric { + template + static void get(MetricType& metric_field, const unsigned long iPoint, + const unsigned short iSensor, MatrixType& mat, unsigned short nDim) { + switch(nDim) { + case 2: { + mat[0][0] = metric_field(iPoint, 0); mat[0][1] = metric_field(iPoint, 1); + mat[1][0] = metric_field(iPoint, 1); mat[1][1] = metric_field(iPoint, 2); + break; + } + case 3: { + mat[0][0] = metric_field(iPoint, 0); mat[0][1] = metric_field(iPoint, 1); mat[0][2] = metric_field(iPoint, 2); + mat[1][0] = metric_field(iPoint, 1); mat[1][1] = metric_field(iPoint, 3); mat[1][2] = metric_field(iPoint, 4); + mat[2][0] = metric_field(iPoint, 2); mat[2][1] = metric_field(iPoint, 4); mat[2][2] = metric_field(iPoint, 5); + break; + } + } + } + + template + static void set(MetricType& metric_field, const unsigned long iPoint, + const unsigned short iSensor, MatrixType& mat, ScalarType scale, + unsigned short nDim) { + switch(nDim) { + case 2: { + metric_field(iPoint, 0) = mat[0][0] * scale; + metric_field(iPoint, 1) = mat[0][1] * scale; + metric_field(iPoint, 2) = mat[1][1] * scale; + break; + } + case 3: { + metric_field(iPoint, 0) = mat[0][0] * scale; + metric_field(iPoint, 1) = mat[0][1] * scale; + metric_field(iPoint, 2) = mat[0][2] * scale; + metric_field(iPoint, 3) = mat[1][1] * scale; + metric_field(iPoint, 4) = mat[1][2] * scale; + metric_field(iPoint, 5) = mat[2][2] * scale; + break; + } + } + } + }; + + struct hessian { + template + static void get(MetricType& metric_field, const unsigned long iPoint, + const unsigned short iSensor, MatrixType& mat, unsigned short nDim) { + switch(nDim) { + case 2: { + mat[0][0] = metric_field(iPoint, iSensor, 0); mat[0][1] = metric_field(iPoint, iSensor, 1); + mat[1][0] = metric_field(iPoint, iSensor, 1); mat[1][1] = metric_field(iPoint, iSensor, 2); + break; + } + case 3: { + mat[0][0] = metric_field(iPoint, iSensor, 0); mat[0][1] = metric_field(iPoint, iSensor, 1); mat[0][2] = metric_field(iPoint, iSensor, 2); + mat[1][0] = metric_field(iPoint, iSensor, 1); mat[1][1] = metric_field(iPoint, iSensor, 3); mat[1][2] = metric_field(iPoint, iSensor, 4); + mat[2][0] = metric_field(iPoint, iSensor, 2); mat[2][1] = metric_field(iPoint, iSensor, 4); mat[2][2] = metric_field(iPoint, iSensor, 5); + break; + } + } + } + + template + static void set(MetricType& metric_field, const unsigned long iPoint, + const unsigned short iSensor, MatrixType& mat, ScalarType scale, + unsigned short nDim) { + switch(nDim) { + case 2: { + metric_field(iPoint, iSensor, 0) = mat[0][0] * scale; + metric_field(iPoint, iSensor, 1) = mat[0][1] * scale; + metric_field(iPoint, iSensor, 2) = mat[1][1] * scale; + break; + } + case 3: { + metric_field(iPoint, iSensor, 0) = mat[0][0] * scale; + metric_field(iPoint, iSensor, 1) = mat[0][1] * scale; + metric_field(iPoint, iSensor, 2) = mat[0][2] * scale; + metric_field(iPoint, iSensor, 3) = mat[1][1] * scale; + metric_field(iPoint, iSensor, 4) = mat[1][2] * scale; + metric_field(iPoint, iSensor, 5) = mat[2][2] * scale; + break; + } + } + } + }; +} + +namespace detail { + +/*! + * \brief Compute determinant of eigenvalues for different dimensions. + * \param[in] EigVal - Array of eigenvalues. + * \return Determinant value. + */ +template +ScalarType computeDeterminant(const ScalarType* EigVal) { + if constexpr (nDim == 2) { + return EigVal[0] * EigVal[1]; + } else if constexpr (nDim == 3) { + return EigVal[0] * EigVal[1] * EigVal[2]; + } else { + static_assert(nDim == 2 || nDim == 3, "Only 2D and 3D supported"); + return ScalarType(0.0); + } +} + +/*! + * \brief Make the eigenvalues of the metrics positive. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in,out] metric - Metric container. + */ +template +void setPositiveDefiniteMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, MetricType& metric) { + + const unsigned long nPointDomain = geometry.GetnPointDomain(); + + ScalarType A[nDim][nDim], EigVec[nDim][nDim], EigVal[nDim], work[nDim]; + + /*--- Minimum eigenvalue threshold ---*/ + const ScalarType eps = 1e-12; + + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + /*--- Get full metric tensor ---*/ + Tensor::get(metric, iPoint, iSensor, A, nDim); + + /*--- Compute eigenvalues and eigenvectors ---*/ + CBlasStructure::EigenDecomposition(A, EigVec, EigVal, nDim, work); + + /*--- Make positive definite by taking absolute value of eigenvalues ---*/ + /*--- Handle NaN and very small values that could cause numerical issues ---*/ + for (auto iDim = 0; iDim < nDim; iDim++) { + if (EigVal[iDim] != EigVal[iDim]) { + /*--- NaN detected, set to small positive value ---*/ + EigVal[iDim] = eps; + } else { + /*--- Take absolute value and ensure minimum threshold ---*/ + EigVal[iDim] = max(fabs(EigVal[iDim]), eps); + } + } + + CBlasStructure::EigenRecomposition(A, EigVec, EigVal, nDim); + + /*--- Store upper half of metric tensor ---*/ + Tensor::set(metric, iPoint, iSensor, A, 1.0, nDim); + } +} + +/*! + * \brief Integrate the Hessian field for the Lp-norm normalization of the metric. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in] metric - Metric container. + * \return Integral of the metric tensor determinant. +*/ +template +ScalarType integrateMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, MetricType& metric) { + + const unsigned long nPointDomain = geometry.GetnPointDomain(); + + /*--- Constants defining normalization ---*/ + const ScalarType p = config.GetMetric_Norm(); + const ScalarType normExp = p / (2.0 * p + nDim); + + ScalarType localIntegral = 0.0; + ScalarType globalIntegral = 0.0; + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + auto nodes = geometry.nodes; + + /*--- Calculate determinant ---*/ + ScalarType det; + if constexpr (nDim == 2) { + const ScalarType m00 = metric(iPoint, 0); + const ScalarType m01 = metric(iPoint, 1); + const ScalarType m11 = metric(iPoint, 2); + det = m00 * m11 - m01 * m01; + } else if constexpr (nDim == 3) { + const ScalarType m00 = metric(iPoint, 0); + const ScalarType m01 = metric(iPoint, 1); + const ScalarType m02 = metric(iPoint, 2); + const ScalarType m11 = metric(iPoint, 3); + const ScalarType m12 = metric(iPoint, 4); + const ScalarType m22 = metric(iPoint, 5); + det = m00 * (m11 * m22 - m12 * m12) - m01 * (m01 * m22 - m02 * m12) + m02 * (m01 * m12 - m02 * m11); + } + + /*--- Integrate determinant ---*/ + const ScalarType Vol = SU2_TYPE::GetValue(nodes->GetVolume(iPoint)); + localIntegral += pow(abs(det), normExp) * Vol; + } + + CBaseMPIWrapper::Allreduce(&localIntegral, &globalIntegral, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + return globalIntegral; +} + +/*! + * \brief Perform an Lp-norm normalization of the metric. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in] integral - Integral of the metric tensor determinant. + * \param[in,out] metric - Metric container. + */ +template +void normalizeMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, ScalarType integral, + MetricType& metric) { + + const unsigned long nPointDomain = geometry.GetnPointDomain(); + + const bool goal = (config.GetGoal_Oriented_Metric()); + + /*--- Constants defining normalization ---*/ + const ScalarType p = config.GetMetric_Norm(); + const ScalarType N = SU2_TYPE::GetValue(config.GetMetric_Complexity()); + const ScalarType globalFactor = pow(N / integral, 2.0 / nDim); + const ScalarType normExp = -1.0 / (2.0 * p + nDim); + + /*--- Size constraints ---*/ + const ScalarType hmin = SU2_TYPE::GetValue(config.GetMetric_Hmin()); + const ScalarType hmax = SU2_TYPE::GetValue(config.GetMetric_Hmax()); + const ScalarType eigmax = 1.0 / pow(hmin, 2.0); + const ScalarType eigmin = 1.0 / pow(hmax, 2.0); + const ScalarType armax2 = pow(SU2_TYPE::GetValue(config.GetMetric_ARmax()), 2.0); + + ScalarType A[nDim][nDim], EigVec[nDim][nDim], EigVal[nDim], work[nDim]; + + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + auto nodes = geometry.nodes; + + /*--- Decompose metric ---*/ + Tensor::get(metric, iPoint, iSensor, A, nDim); + CBlasStructure::EigenDecomposition(A, EigVec, EigVal, nDim, work); + + /*--- Normalize eigenvalues ---*/ + const ScalarType det = computeDeterminant(EigVal); + const ScalarType factor = globalFactor * pow(abs(det), normExp); + for (auto iDim = 0u; iDim < nDim; ++iDim) + EigVal[iDim] = factor * EigVal[iDim]; + + /*--- Clip by user-specified size constraints ---*/ + for (auto iDim = 0u; iDim < nDim; ++iDim) + EigVal[iDim] = min(max(abs(EigVal[iDim]), eigmin), eigmax); + + /*--- Clip by user-specified aspect ratio ---*/ + unsigned short iMax = 0; + for (auto iDim = 1; iDim < nDim; ++iDim) + iMax = (EigVal[iDim] > EigVal[iMax])? iDim : iMax; + + for (auto iDim = 0u; iDim < nDim; ++iDim) + EigVal[iDim] = max(EigVal[iDim], EigVal[iMax]/armax2); + + /*--- Recompose and store metric ---*/ + CBlasStructure::EigenRecomposition(A, EigVec, EigVal, nDim); + Tensor::set(metric, iPoint, iSensor, A, 1.0, nDim); + } +} +} // end namespace detail + +/*! + * \brief Make the eigenvalues of the metrics positive. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in,out] metric - Metric container. + */ +template +void setPositiveDefiniteMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, MetricType& metric) { + switch (geometry.GetnDim()) { + case 2: + detail::setPositiveDefiniteMetrics<2, ScalarType, Tensor>(geometry, config, iSensor, metric); + break; + case 3: + detail::setPositiveDefiniteMetrics<3, ScalarType, Tensor>(geometry, config, iSensor, metric); + break; + default: + SU2_MPI::Error("Too many dimensions for metric computation.", CURRENT_FUNCTION); + break; + } +} + +/*! + * \brief Integrate the Hessian field for the Lp-norm normalization of the metric. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in] metric - Metric container. + * \return Integral of the metric tensor determinant. +*/ +template +ScalarType integrateMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, MetricType& metric) { + su2double integral; + switch (geometry.GetnDim()) { + case 2: + integral = detail::integrateMetrics<2, ScalarType>(geometry, config, iSensor, metric); + break; + case 3: + integral = detail::integrateMetrics<3, ScalarType>(geometry, config, iSensor, metric); + break; + default: + SU2_MPI::Error("Too many dimensions for metric integration.", CURRENT_FUNCTION); + break; + } + + return integral; +} + +/*! + * \brief Perform an Lp-norm normalization of the metric. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + * \param[in] integral - Integral of the metric tensor determinant. + * \param[in,out] metric - Metric container. + */ +template +void normalizeMetrics(CGeometry& geometry, const CConfig& config, + unsigned short iSensor, ScalarType integral, + MetricType& metric) { + switch (geometry.GetnDim()) { + case 2: + detail::normalizeMetrics<2, ScalarType, Tensor>(geometry, config, iSensor, integral, metric); + break; + case 3: + detail::normalizeMetrics<3, ScalarType, Tensor>(geometry, config, iSensor, integral, metric); + break; + default: + SU2_MPI::Error("Too many dimensions for metric normalization.", CURRENT_FUNCTION); + break; + } +} diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index e3f029a37538..1883e8e26d2a 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -273,6 +273,49 @@ class CEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::MACH: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetVelocity2(iPoint) / nodes->GetSoundSpeed(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::DENSITY: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetDensity(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + } + } + } + else { + // TODO: store variables which need Hessian computation for goal-oriented metric + } + } + /*! * \brief Set gradients of coefficients for fixed CL mode * \param[in] config - Definition of the particular problem. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index f0b3a4c493d1..9f801d22925a 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -85,7 +85,8 @@ class CSolver { nSecondaryVar, /*!< \brief Number of primitive variables of the problem. */ nSecondaryVarGrad, /*!< \brief Number of primitive variables of the problem in the gradient computation. */ nVarGrad, /*!< \brief Number of variables for deallocating the LS Cvector. */ - nDim; /*!< \brief Number of dimensions of the problem. */ + nDim, /*!< \brief Number of dimensions of the problem. */ + nSymMat; /*!< \brief Number of symmetric matrix componenents for Hessian and metric tensor. */ unsigned long nPoint; /*!< \brief Number of points of the computational grid. */ unsigned long nPointDomain; /*!< \brief Number of points of the computational grid. */ su2double Max_Delta_Time, /*!< \brief Maximum value of the delta time for all the control volumes. */ @@ -580,6 +581,22 @@ class CSolver { */ inline virtual void SetPrimitive_Limiter(CGeometry *geometry, const CConfig *config) { } + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + virtual void SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) { } + + /*! + * \brief Compute the Green-Gauss Hessian of the solution. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] idxVel - Index to velocity, -1 if no velocity is present in the solver. + * \param[in] reconstruction - indicator that the gradient being computed is for upwind reconstruction. + */ + void SetHessian_GG(CGeometry *geometry, const CConfig *config, short idxVel, const unsigned short Kind_Solver); + /*! * \brief Compute the projection of a variable for MUSCL reconstruction. * \note The result should be halved when added to i (or subtracted from j). @@ -4349,24 +4366,42 @@ class CSolver { END_SU2_OMP_FOR } -inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container, - CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { + inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container, + CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { - AD::StartNoSharedReading(); + AD::StartNoSharedReading(); - SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads())) - for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { - /*--- Get control volume size. ---*/ - su2double Volume = geometry->nodes->GetVolume(iPoint); - /*--- Compute the residual for this control volume and subtract. ---*/ - for (auto iVar = 0ul; iVar < nVar; iVar++) { - LinSysRes(iPoint,iVar) -= base_nodes->GetUserDefinedSource()(iPoint, iVar) * Volume; + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + /*--- Get control volume size. ---*/ + su2double Volume = geometry->nodes->GetVolume(iPoint); + /*--- Compute the residual for this control volume and subtract. ---*/ + for (auto iVar = 0ul; iVar < nVar; iVar++) { + LinSysRes(iPoint,iVar) -= base_nodes->GetUserDefinedSource()(iPoint, iVar) * Volume; + } } + END_SU2_OMP_FOR + + AD::EndNoSharedReading(); } - END_SU2_OMP_FOR - AD::EndNoSharedReading(); -} + /*! + * \brief Compute the metric tensor field. + * \param[in] solver - Physical definition of the problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig *config); + + /*! + * \brief Add contribution to the metric field. + * \param[in] solver - Physical definition of the problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iSensor - Index of the sensor to work on. + */ + void AddMetrics(CSolver **solver, const CGeometry *geometry, const CConfig *config, + const unsigned short iSensor); protected: /*! diff --git a/SU2_CFD/src/solvers/CBaselineSolver.cpp b/SU2_CFD/src/solvers/CBaselineSolver.cpp index fa89ee5e4830..ef6703f7efb2 100644 --- a/SU2_CFD/src/solvers/CBaselineSolver.cpp +++ b/SU2_CFD/src/solvers/CBaselineSolver.cpp @@ -40,6 +40,7 @@ CBaselineSolver::CBaselineSolver(CGeometry *geometry, CConfig *config) { /*--- Define geometry constants in the solver structure ---*/ nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); /*--- Routines to access the number of variables and string names. ---*/ diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 8bae5e6ca36f..2e8a1b3a0aaa 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -41,6 +41,7 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di nVar = direct_solver->GetnVar(); nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); nMarker = config->GetnMarker_All(); nPoint = geometry->GetnPoint(); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index dcca14edd19d..ebb02534b492 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -114,6 +114,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, ---*/ nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); nVar = nDim + 2; nPrimVar = nDim + 9; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 3df1ffbe3dd8..53d672829298 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -104,6 +104,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned * Incompressible flow, primitive variables (P, vx, vy, vz, T, rho, beta, lamMu, EddyMu, Kt_eff, Cp, Cv) ---*/ nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); /*--- Make sure to align the sizes with the constructor of CIncEulerVariable. ---*/ nVar = nDim + 2; diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 65f3ea81f123..1497f4fe0f70 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -92,6 +92,7 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config, nSpecies = config->GetnSpecies(); nMarker = config->GetnMarker_All(); nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); nPoint = geometry->GetnPoint(); nPointDomain = geometry->GetnPointDomain(); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index a9c3875406bf..16bbfc1ea3b1 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -29,6 +29,7 @@ #include "../../include/solvers/CSolver.hpp" #include "../../include/gradients/computeGradientsGreenGauss.hpp" #include "../../include/gradients/computeGradientsLeastSquares.hpp" +#include "../../include/gradients/computeMetrics.hpp" #include "../../include/limiters/computeLimiters.hpp" #include "../../../Common/include/toolboxes/MMS/CIncTGVSolution.hpp" #include "../../../Common/include/toolboxes/MMS/CInviscidVortexSolution.hpp" @@ -1379,6 +1380,14 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; + case MPI_QUANTITIES::GRADIENT_ADAPT: + COUNT_PER_POINT = config->GetGoal_Oriented_Metric()? nVar*nDim : config->GetnMetric_Sensor()*nDim; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; + case MPI_QUANTITIES::HESSIAN: + COUNT_PER_POINT = config->GetGoal_Oriented_Metric()? nVar*nSymMat : config->GetnMetric_Sensor()*nSymMat; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; default: SU2_MPI::Error("Unrecognized quantity for point-to-point MPI comms.", CURRENT_FUNCTION); @@ -1393,6 +1402,8 @@ namespace CommHelpers { case MPI_QUANTITIES::PRIMITIVE_GRADIENT: return nodes->GetGradient_Primitive(); case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: return nodes->GetGradient_Reconstruction(); case MPI_QUANTITIES::AUXVAR_GRADIENT: return nodes->GetAuxVarGradient(); + case MPI_QUANTITIES::GRADIENT_ADAPT: return nodes->GetGradient_Adapt(); + case MPI_QUANTITIES::HESSIAN: return nodes->GetHessian(); default: return nodes->GetGradient(); } } @@ -1409,7 +1420,7 @@ void CSolver::InitiateComms(CGeometry *geometry, /*--- Local variables ---*/ - unsigned short iVar, iDim; + unsigned short iVar, iDim, iMat; unsigned short COUNT_PER_POINT = 0; unsigned short MPI_TYPE = 0; @@ -1435,6 +1446,7 @@ void CSolver::InitiateComms(CGeometry *geometry, /*--- Handle the different types of gradient and limiter. ---*/ const auto nVarGrad = COUNT_PER_POINT / nDim; + const auto nVarHess = COUNT_PER_POINT / nSymMat; auto& gradient = CommHelpers::selectGradient(base_nodes, commType); auto& limiter = CommHelpers::selectLimiter(base_nodes, commType); @@ -1503,10 +1515,16 @@ void CSolver::InitiateComms(CGeometry *geometry, case MPI_QUANTITIES::SOLUTION_GRAD_REC: case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: case MPI_QUANTITIES::AUXVAR_GRADIENT: + case MPI_QUANTITIES::GRADIENT_ADAPT: for (iVar = 0; iVar < nVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset+iVar*nDim+iDim] = gradient(iPoint, iVar, iDim); break; + case MPI_QUANTITIES::HESSIAN: + for (iVar = 0; iVar < nVarHess; iVar++) + for (iMat = 0; iMat < nSymMat; iMat++) + bufDSend[buf_offset+iVar*nSymMat+iMat] = gradient(iPoint, iVar, iMat); + break; case MPI_QUANTITIES::SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); @@ -1551,7 +1569,7 @@ void CSolver::CompleteComms(CGeometry *geometry, /*--- Local variables ---*/ - unsigned short iDim, iVar; + unsigned short iDim, iVar, iMat; unsigned long iPoint, iRecv, nRecv, msg_offset, buf_offset; unsigned short COUNT_PER_POINT = 0; unsigned short MPI_TYPE = 0; @@ -1572,6 +1590,7 @@ void CSolver::CompleteComms(CGeometry *geometry, /*--- Handle the different types of gradient and limiter. ---*/ const auto nVarGrad = COUNT_PER_POINT / nDim; + const auto nVarHess = COUNT_PER_POINT / nSymMat; auto& gradient = CommHelpers::selectGradient(base_nodes, commType); auto& limiter = CommHelpers::selectLimiter(base_nodes, commType); @@ -1651,10 +1670,16 @@ void CSolver::CompleteComms(CGeometry *geometry, case MPI_QUANTITIES::SOLUTION_GRAD_REC: case MPI_QUANTITIES::PRIMITIVE_GRAD_REC: case MPI_QUANTITIES::AUXVAR_GRADIENT: + case MPI_QUANTITIES::GRADIENT_ADAPT: for (iVar = 0; iVar < nVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) gradient(iPoint,iVar,iDim) = bufDRecv[buf_offset+iVar*nDim+iDim]; break; + case MPI_QUANTITIES::HESSIAN: + for (iVar = 0; iVar < nVarHess; iVar++) + for (iMat = 0; iMat < nSymMat; iMat++) + gradient(iPoint, iVar, iMat) = bufDRecv[buf_offset+iVar*nSymMat+iMat]; + break; case MPI_QUANTITIES::SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); @@ -2165,6 +2190,20 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, const CConfig *config computeGradientsLeastSquares(this, comm, commPer, *geometry, *config, weighted, solution, 0, nVar, idxVel, gradient, rmatrix); } +void CSolver::SetHessian_GG(CGeometry *geometry, const CConfig *config, short idxVel, const unsigned short Kind_Solver) { + const auto& solution = config->GetGoal_Oriented_Metric()? base_nodes->GetSolution() : base_nodes->GetPrimitive_Adapt(); + auto& gradient = base_nodes->GetGradient_Adapt(); + auto nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + + computeGradientsGreenGauss(this, MPI_QUANTITIES::GRADIENT_ADAPT, PERIODIC_GRAD_ADAPT, + *geometry, *config, solution, 0, nHess, idxVel, gradient); + + auto& hessian = base_nodes->GetHessian(); + + computeHessiansGreenGauss(this, MPI_QUANTITIES::HESSIAN, PERIODIC_HESSIAN, + *geometry, *config, gradient, 0, nHess, idxVel, hessian); +} + void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) { /*--- Loop domain points. ---*/ @@ -4346,3 +4385,93 @@ void CSolver::SavelibROM(CGeometry *geometry, CConfig *config, bool converged) { #endif } + + +void CSolver::ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig *config) { + /*--- TODO: - goal-oriented metric ---*/ + /*--- - metric intersection ---*/ + const unsigned long nPointDomain = geometry->GetnPointDomain(); + + const bool visc = (config->GetViscous()); + const bool turb = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); + + const bool goal = (config->GetGoal_Oriented_Metric()); + const bool normalize = (config->GetNormalize_Metric()); + + unsigned short nSensor = config->GetnMetric_Sensor(); + + const unsigned long time_iter = config->GetTimeIter(); + const bool steady = (config->GetTime_Marching() == TIME_MARCHING::STEADY); + const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) || + (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); + const bool is_last_iter = (time_iter == config->GetnTime_Iter() - 1) || (steady); + + /*--- Integrate and normalize the metric tensor field ---*/ + vector integrals; + for (auto iSensor = 0u; iSensor < nSensor; ++iSensor) { + SU2_OMP_MASTER + + if (!goal) { + /*--- Make the Hessian eigenvalues positive definite, and add to the metric tensor ---*/ + auto& hessians = base_nodes->GetHessian(); + setPositiveDefiniteMetrics(*geometry, *config, iSensor, hessians); + AddMetrics(solver, geometry, config, iSensor); + + /*--- Integrate metric field on the last iteration (the end of the simulation if steady) ---*/ + auto& metrics = base_nodes->GetMetric(); + double integral = 0.0; + if (is_last_iter) + integral = integrateMetrics(*geometry, *config, iSensor, metrics); + + /*--- Normalize the metric field for steady simulations, or if requested for unsteady ---*/ + if (steady || (normalize && is_last_iter)) + normalizeMetrics(*geometry, *config, iSensor, integral, metrics); + + /*--- Store the integral to be written ---*/ + if (is_last_iter) { + integrals.push_back(integral); + if (rank == MASTER_NODE) { + cout << "Global metric normalization integral for sensor "; + cout << config->GetMetric_SensorString(iSensor) << ": " << integral << endl; + } + } + } else { + SU2_MPI::Error("Goal-oriented metric not currently implemented.", CURRENT_FUNCTION); + } + END_SU2_OMP_MASTER + SU2_OMP_BARRIER + } +} + +void CSolver::AddMetrics(CSolver **solver, const CGeometry*geometry, const CConfig *config, + const unsigned short iSensor) { + /*--- TODO: - goal-oriented metric ---*/ + /*--- - metric intersection ---*/ + auto varFlo = solver[FLOW_SOL]->GetNodes(); + + const unsigned long nPointDomain = geometry->GetnPointDomain(); + const unsigned short nSymMat = 3*(nDim-1); + const unsigned short nVarFlo = solver[FLOW_SOL]->GetnVar(); + const unsigned short nSensor = config->GetnMetric_Sensor(); + + const bool turb = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); + const bool goal = (config->GetGoal_Oriented_Metric()); + + const unsigned long time_iter = config->GetTimeIter(); + const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) || + (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); + const bool is_first_iter = (time_iter == 0); + const bool is_last_iter = (time_iter == config->GetnTime_Iter() - 1); + + double coeff = (time_stepping && (is_first_iter || is_last_iter))? 0.5 : 1.0; + if (time_stepping) coeff *= SU2_TYPE::GetValue(config->GetTime_Step()); + + for(auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + for (auto iMat = 0; iMat < nSymMat; ++iMat) { + double hess = SU2_TYPE::GetValue(varFlo->GetHessian(iPoint, iSensor, iMat)); + varFlo->AddMetric(iPoint, iMat, coeff * hess); + } + } +} diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 7649165d0be3..0124d6e887b5 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -53,6 +53,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor /*--- Define geometry constants in the solver structure ---*/ nDim = geometry->GetnDim(); + nSymMat = 3 * (nDim - 1); /*--- Single grid simulation ---*/ @@ -301,7 +302,7 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry* /*--- Roughness heights. ---*/ numerics->SetRoughness(geometry->nodes->GetRoughnessHeight(iPoint), geometry->nodes->GetRoughnessHeight(jPoint)); }; - + /*--- Now instantiate the generic non-conservative implementation with the functor above. ---*/ Viscous_Residual_NonCons(iEdge, geometry, solver_container, numerics, config, SolverSpecificNumerics); From 14aba45873b5d1e2819a4f4714c4ffe0ae319d94 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 12:03:16 -0800 Subject: [PATCH 04/12] Add single zone driver calls to calculate metric fiedl --- SU2_CFD/include/drivers/CSinglezoneDriver.hpp | 5 +++ SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 35 +++++++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp index 355369556ea5..d7e789ae66ce 100644 --- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp +++ b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp @@ -110,4 +110,9 @@ class CSinglezoneDriver : public CDriver { */ bool Monitor(unsigned long TimeIter) override; + /*! + * \brief Perform all steps to compute the metric tensor. + */ + virtual void ComputeMetricField(); + }; diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 1270eef237ab..eec6ec19bc9a 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -175,6 +175,11 @@ void CSinglezoneDriver::Postprocess() { iteration_container[ZONE_0][INST_0]->Relaxation(output_container[ZONE_0], integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + /*--- Compute metric for anisotropic mesh adaptation ---*/ + + if (config_container[ZONE_0]->GetCompute_Metric()) + ComputeMetricField(); + } void CSinglezoneDriver::Update() { @@ -312,3 +317,33 @@ bool CSinglezoneDriver::Monitor(unsigned long TimeIter){ bool CSinglezoneDriver::GetTimeConvergence() const{ return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]); } + +void CSinglezoneDriver::ComputeMetricField() { + + auto solver = solver_container[ZONE_0][INST_0][MESH_0]; + auto solver_flow = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; + auto geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + auto config = config_container[ZONE_0]; + int idxVel = -1; + + if (rank == MASTER_NODE){ + cout << endl <<"----------------------------- Compute Metric ----------------------------" << endl; + cout << "Storing primitive variables needed for gradients in metric." << endl; + } + solver_flow->InitiateComms(geometry, config, MPI_QUANTITIES::SOLUTION); + solver_flow->CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); + solver_flow->Preprocessing(geometry, solver, config, MESH_0, NO_RK_ITER, + RUNTIME_FLOW_SYS, true); + solver_flow->SetPrimitive_Adapt(geometry, config, solver_flow->GetNodes()); + + if (config->GetKind_Hessian_Method() == GREEN_GAUSS) { + if(rank == MASTER_NODE) cout << "Computing Hessians using Green-Gauss." << endl; + solver_flow->SetHessian_GG(geometry, config, idxVel, RUNTIME_FLOW_SYS); + } + else { + SU2_MPI::Error("Unsupported Hessian method.", CURRENT_FUNCTION); + } + + if(rank == MASTER_NODE) cout << "Computing feature-based metric tensor." << endl; + solver_flow->ComputeMetric(solver, geometry, config); +} From 631f4b57bf6c98305786f9fbac87a94ee095f0bc Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 12:10:26 -0800 Subject: [PATCH 05/12] Add metric tensor outputs --- SU2_CFD/include/output/CFlowOutput.hpp | 16 +++++++++++ SU2_CFD/src/output/CFlowCompOutput.cpp | 3 ++ SU2_CFD/src/output/CFlowOutput.cpp | 39 ++++++++++++++++++++++++++ 3 files changed, 58 insertions(+) diff --git a/SU2_CFD/include/output/CFlowOutput.hpp b/SU2_CFD/include/output/CFlowOutput.hpp index 63867b091ae4..7c8c23af181f 100644 --- a/SU2_CFD/include/output/CFlowOutput.hpp +++ b/SU2_CFD/include/output/CFlowOutput.hpp @@ -352,4 +352,20 @@ class CFlowOutput : public CFVMOutput{ * \param[in] config - Definition of the particular problem per zone. */ void SetFixedCLScreenOutput(const CConfig *config); + + /*! + * \brief Add mesh adaptation outputs. + * \param[in] config - Definition of the particular problem. + */ + void AddMeshAdaptationOutputs(const CConfig* config); + + /*! + * \brief Load mesh adaptation outputs. + * \param[in] config - Definition of the particular problem. + * \param[in] solver - The container holding all solution data. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] iPoint - Index of the point. + */ + void LoadMeshAdaptationOutputs(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry, + unsigned long iPoint); }; diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 2476498e84d7..6bcae5a287f0 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -303,6 +303,9 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ AddCommonFVMOutputs(config); + // Metric tensor and dual-cell volume + AddMeshAdaptationOutputs(config); + if (config->GetTime_Domain()) { SetTimeAveragedFields(); } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index b355c0698afb..be5647734a0b 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -4109,3 +4109,42 @@ void CFlowOutput::AddTurboOutput(unsigned short nZone){ AddHistoryOutput("KineticEnergyLoss_Stage", "KELC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Kinetic Energy Loss Coefficient", HistoryFieldType::DEFAULT); AddHistoryOutput("TotPressureLoss_Stage", "TPLC_all", ScreenOutputFormat::SCIENTIFIC, "TURBO_PERF", "Machine Pressure Loss Coefficient", HistoryFieldType::DEFAULT); } + +void CFlowOutput::AddMeshAdaptationOutputs(const CConfig* config) { + + // Anisotropic metric tensor, and dual-cell volume + if(config->GetCompute_Metric()) { + // Common metric components for both 2D and 3D + AddVolumeOutput("METRIC_XX", "Metric_xx", "MESH_ADAPT", "x-x-component of the metric"); + AddVolumeOutput("METRIC_XY", "Metric_xy", "MESH_ADAPT", "x-y-component of the metric"); + AddVolumeOutput("METRIC_YY", "Metric_yy", "MESH_ADAPT", "y-y-component of the metric"); + + // Additional components for 3D + if (nDim == 3) { + AddVolumeOutput("METRIC_XZ", "Metric_xz", "MESH_ADAPT", "x-z-component of the metric"); + AddVolumeOutput("METRIC_YZ", "Metric_yz", "MESH_ADAPT", "y-z-component of the metric"); + AddVolumeOutput("METRIC_ZZ", "Metric_zz", "MESH_ADAPT", "z-z-component of the metric"); + } + } +} + +void CFlowOutput::LoadMeshAdaptationOutputs(const CConfig* config, const CSolver* const* solver, const CGeometry* geometry, + unsigned long iPoint) { + + const auto* Node_Flow = solver[FLOW_SOL]->GetNodes(); + const auto* Node_Geo = geometry->nodes; + if(config->GetCompute_Metric()) { + // Common metric components for both 2D and 3D + SetVolumeOutputValue("METRIC_XX", iPoint, Node_Flow->GetMetric(iPoint, 0)); + SetVolumeOutputValue("METRIC_XY", iPoint, Node_Flow->GetMetric(iPoint, 1)); + SetVolumeOutputValue("METRIC_YY", iPoint, Node_Flow->GetMetric(iPoint, 2)); + + // Additional components for 3D + if (nDim == 3) { + SetVolumeOutputValue("METRIC_XZ", iPoint, Node_Flow->GetMetric(iPoint, 3)); + SetVolumeOutputValue("METRIC_YZ", iPoint, Node_Flow->GetMetric(iPoint, 4)); + SetVolumeOutputValue("METRIC_ZZ", iPoint, Node_Flow->GetMetric(iPoint, 5)); + } + } +} + From 14dddf6211ed5c6fdd9c61f6e640a9acafb78ce9 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 12:11:56 -0800 Subject: [PATCH 06/12] Add tensor outputs to paraview filewriter --- .../filewriter/CParaviewXMLFileWriter.cpp | 178 ++++++++++++++++-- 1 file changed, 158 insertions(+), 20 deletions(-) diff --git a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp index 5a848b2e11a9..14fd3e1a3fa2 100644 --- a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp +++ b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp @@ -53,8 +53,10 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){ } /*--- We always have 3 coords, independent of the actual value of nDim ---*/ + /*--- We also always have 6 tensor components ---*/ const int NCOORDS = 3; + const int NTENSOR = 6; const unsigned short nDim = dataSorter->GetnDim(); unsigned short iDim = 0; @@ -139,30 +141,85 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){ fieldname.erase(remove(fieldname.begin(), fieldname.end(), '"'), fieldname.end()); - /*--- Check whether this field is a vector or scalar. ---*/ + /*--- Check whether this field is a vector, tensor, or scalar. ---*/ + /*--- Check tensor components first (longer patterns) to avoid false matches ---*/ - bool output_variable = true, isVector = false; - size_t found = fieldNames[iField].find("_x"); + bool output_variable = true, isVector = false, isTensor = false; + + // Check for tensor components first (_xx, _xy, _yy, _xz, _yz, _zz) + size_t found = fieldNames[iField].find("_xx"); if (found!=string::npos) { output_variable = true; - isVector = true; + isTensor = true; + } + found = fieldNames[iField].find("_xy"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the XY component. ---*/ + output_variable = false; + isTensor = true; + VarCounter++; + } + found = fieldNames[iField].find("_yy"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the YY component. ---*/ + output_variable = false; + isTensor = true; + VarCounter++; + } + found = fieldNames[iField].find("_xz"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the XZ component. ---*/ + output_variable = false; + isTensor = true; + VarCounter++; } - found = fieldNames[iField].find("_y"); + found = fieldNames[iField].find("_yz"); if (found!=string::npos) { - /*--- We have found a vector, so skip the Y component. ---*/ + /*--- We have found a tensor, so skip the YZ component. ---*/ output_variable = false; + isTensor = true; VarCounter++; } - found = fieldNames[iField].find("_z"); + found = fieldNames[iField].find("_zz"); if (found!=string::npos) { - /*--- We have found a vector, so skip the Z component. ---*/ + /*--- We have found a tensor, so skip the ZZ component. ---*/ output_variable = false; + isTensor = true; VarCounter++; } - /*--- Write the point data as an vector or a scalar. ---*/ + // Check for vector components only if not a tensor (_x, _y, _z) + if (!isTensor) { + found = fieldNames[iField].find("_x"); + if (found!=string::npos) { + output_variable = true; + isVector = true; + } + found = fieldNames[iField].find("_y"); + if (found!=string::npos) { + /*--- We have found a vector, so skip the Y component. ---*/ + output_variable = false; + VarCounter++; + } + found = fieldNames[iField].find("_z"); + if (found!=string::npos) { + /*--- We have found a vector, so skip the Z component. ---*/ + output_variable = false; + VarCounter++; + } + } + + /*--- Write the point data as an vector, tensor, or a scalar. ---*/ + + if (output_variable && isTensor) { + + /*--- Adjust the string name to remove the tensor component suffix ---*/ - if (output_variable && isVector) { + fieldname.erase(fieldname.end()-3,fieldname.end()); + + AddDataArray(VTKDatatype::FLOAT32, fieldname, NTENSOR, myPoint*NTENSOR, GlobalPoint*NTENSOR); + + } else if (output_variable && isVector) { /*--- Adjust the string name to remove the leading "X-" ---*/ @@ -252,33 +309,112 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){ VarCounter = varStart; for (iField = varStart; iField < fieldNames.size(); iField++) { - /*--- Check whether this field is a vector or scalar. ---*/ + /*--- Check whether this field is a vector, tensor, or scalar. ---*/ + /*--- Check tensor components first (longer patterns) to avoid double counting ---*/ + + bool output_variable = true, isVector = false, isTensor = false; - bool output_variable = true, isVector = false; - size_t found = fieldNames[iField].find("_x"); + // Check for tensor components first (_xx, _xy, _yy, _xz, _yz, _zz) + size_t found = fieldNames[iField].find("_xx"); if (found!=string::npos) { output_variable = true; - isVector = true; + isTensor = true; } - found = fieldNames[iField].find("_y"); + found = fieldNames[iField].find("_xy"); if (found!=string::npos) { - /*--- We have found a vector, so skip the Y component. ---*/ + /*--- We have found a tensor, so skip the XY component. ---*/ output_variable = false; + isTensor = true; VarCounter++; } - found = fieldNames[iField].find("_z"); + found = fieldNames[iField].find("_yy"); if (found!=string::npos) { - /*--- We have found a vector, so skip the Z component. ---*/ + /*--- We have found a tensor, so skip the YY component. ---*/ output_variable = false; + isTensor = true; + VarCounter++; + } + found = fieldNames[iField].find("_xz"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the XZ component. ---*/ + output_variable = false; + isTensor = true; + VarCounter++; + } + found = fieldNames[iField].find("_yz"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the YZ component. ---*/ + output_variable = false; + isTensor = true; + VarCounter++; + } + found = fieldNames[iField].find("_zz"); + if (found!=string::npos) { + /*--- We have found a tensor, so skip the ZZ component. ---*/ + output_variable = false; + isTensor = true; VarCounter++; } - /*--- Write the point data as an vector or a scalar. ---*/ + // Check for vector components only if not a tensor (_x, _y, _z) + if (!isTensor) { + found = fieldNames[iField].find("_x"); + if (found!=string::npos) { + output_variable = true; + isVector = true; + } + found = fieldNames[iField].find("_y"); + if (found!=string::npos) { + /*--- We have found a vector, so skip the Y component. ---*/ + output_variable = false; + VarCounter++; + } + found = fieldNames[iField].find("_z"); + if (found!=string::npos) { + /*--- We have found a vector, so skip the Z component. ---*/ + output_variable = false; + VarCounter++; + } + } - if (output_variable && isVector) { + /*--- Write the point data as a tensor, vector, or scalar. ---*/ + + if (output_variable && isTensor) { + + /*--- Resize buffer to accommodate tensor data ---*/ + dataBufferFloat.resize(myPoint*NTENSOR); + + float val = 0.0; + for (iPoint = 0; iPoint < myPoint; iPoint++) { + dataBufferFloat[iPoint*NTENSOR + 0] = (float)dataSorter->GetData(VarCounter+0,iPoint); // XX + dataBufferFloat[iPoint*NTENSOR + 1] = (float)dataSorter->GetData(VarCounter+2,iPoint); // YY + dataBufferFloat[iPoint*NTENSOR + 3] = (float)dataSorter->GetData(VarCounter+1,iPoint); // XY + if (nDim == 2) { + /*--- For 2D tensors, set the off-diagonal z-components to 0, and zz-component to 1 ---*/ + + dataBufferFloat[iPoint*NTENSOR + 2] = 1.0; // ZZ = 1 + dataBufferFloat[iPoint*NTENSOR + 4] = 0.0; // YZ = 0 + dataBufferFloat[iPoint*NTENSOR + 5] = 0.0; // XZ = 0 + } else { + /*--- For 3D tensors, get the z-components ---*/ + dataBufferFloat[iPoint*NTENSOR + 2] = (float)dataSorter->GetData(VarCounter+5,iPoint); // ZZ + dataBufferFloat[iPoint*NTENSOR + 4] = (float)dataSorter->GetData(VarCounter+4,iPoint); // YZ + dataBufferFloat[iPoint*NTENSOR + 5] = (float)dataSorter->GetData(VarCounter+3,iPoint); // XZ + } + } + + WriteDataArray(dataBufferFloat.data(), VTKDatatype::FLOAT32, myPoint*NTENSOR, GlobalPoint*NTENSOR, + dataSorter->GetnPointCumulative(rank)*NTENSOR); + + VarCounter++; + + } else if (output_variable && isVector) { /*--- Load up the buffer for writing this rank's vector data. ---*/ + /*--- Resize buffer back to coordinate size ---*/ + dataBufferFloat.resize(myPoint*NCOORDS); + float val = 0.0; for (iPoint = 0; iPoint < myPoint; iPoint++) { for (iDim = 0; iDim < NCOORDS; iDim++) { @@ -298,6 +434,8 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){ } else if (output_variable) { + /*--- For scalar data, ensure buffer is properly sized ---*/ + dataBufferFloat.resize(myPoint); /*--- For now, create a temp 1D buffer to load up the data for writing. This will be replaced with a derived data type most likely. ---*/ From c8420370a8e4c52fd513adfdbdf954cc4478cc4f Mon Sep 17 00:00:00 2001 From: bmunguia Date: Mon, 3 Nov 2025 12:32:55 -0800 Subject: [PATCH 07/12] Add Alberto/Andrea's Primitive_Adapt for NEMO and incompressible Euler --- Common/src/CConfig.cpp | 4 +- SU2_CFD/include/solvers/CEulerSolver.hpp | 21 ++++++-- SU2_CFD/include/solvers/CIncEulerSolver.hpp | 49 +++++++++++++++++ SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 55 ++++++++++++++++++++ 4 files changed, 123 insertions(+), 6 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 5232edc31896..a647545b7801 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5814,7 +5814,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } if (Kind_Solver == MAIN_SOLVER::NEMO_EULER || Kind_Solver == MAIN_SOLVER::NEMO_NAVIER_STOKES) { - if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY) + if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY || Metric_Sensor[iSensor] == METRIC_SENSOR::TOTAL_PRESSURE) SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for NEMO problems."), CURRENT_FUNCTION); } if (Kind_Solver != MAIN_SOLVER::NEMO_EULER && Kind_Solver != MAIN_SOLVER::NEMO_NAVIER_STOKES) { @@ -5822,7 +5822,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for non-NEMO problems."), CURRENT_FUNCTION); } if (Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) { - if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::ENERGY) + if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::MACH) SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for INC problems."), CURRENT_FUNCTION); } } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 1883e8e26d2a..0109d461015e 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -283,9 +283,9 @@ class CEulerSolver : public CFVMFlowSolverBaseGetnMetric_Sensor(); for (auto iSensor = 0; iSensor < nSensor; iSensor++) { switch (config->GetMetric_Sensor(iSensor)) { - case METRIC_SENSOR::MACH: + case METRIC_SENSOR::DENSITY: for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetVelocity2(iPoint) / nodes->GetSoundSpeed(iPoint); + const su2double prim_var = nodes->GetDensity(iPoint); nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); } break; @@ -295,16 +295,29 @@ class CEulerSolver : public CFVMFlowSolverBaseSetPrimitive_Adapt(iPoint, iSensor, prim_var); } break; + case METRIC_SENSOR::TOTAL_PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double Mach = sqrt(nodes->GetVelocity2(iPoint))/nodes->GetSoundSpeed(iPoint); + const su2double prim_var = nodes->GetPressure(iPoint)*pow((1+0.2*Mach*Mach), 3.5); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; case METRIC_SENSOR::TEMPERATURE: for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { const su2double prim_var = nodes->GetTemperature(iPoint); nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); } break; - case METRIC_SENSOR::DENSITY: + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::MACH: default: for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetDensity(iPoint); + const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); } break; diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 54fd616d0de8..c9fe7191d6ff 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -86,6 +86,55 @@ class CIncEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::DENSITY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetDensity(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TOTAL_PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint) + 0.5*nodes->GetDensity(iPoint)*nodes->GetVelocity2(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::PRESSURE: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + } + } + } + else { + // TODO: store variables which need Hessian computation for goal-oriented metric + } + } + /*! * \brief Update the Beta parameter for the incompressible preconditioner. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index 9077ea195d89..1f13399554c3 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -195,6 +195,61 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE_VE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature_ve(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY_VE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = *nodes->GetEve(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::MACH: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + } + } + } + else { + // TODO: store variables which need Hessian computation for goal-oriented metric + } + } + /*! * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over * a nonlinear iteration for stability. From 10adbacfbbb26e00dba6378200f9a0f3e6ef9719 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Tue, 4 Nov 2025 10:16:25 -0800 Subject: [PATCH 08/12] Update adaptation python option descriptions --- Common/src/CConfig.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a647545b7801..95a993c46eb6 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3072,15 +3072,15 @@ void CConfig::SetConfig_Options() { addPythonOption("ADAP_ADJ_CFLS"); /*!\brief ADAP_RESIDUAL_REDUCTIONS \n DESCRIPTION: Residual reduction at each target complexity */ addPythonOption("ADAP_RESIDUAL_REDUCTIONS"); - /*!\brief ADAP_ITER \n DESCRIPTION: Maximum cell size at each target complexity */ + /*!\brief ADAP_HMAXS \n DESCRIPTION: Maximum cell size at each target complexity */ addPythonOption("ADAP_HMAXS"); - /*!\brief ADAP_ITER \n DESCRIPTION: Minimum cell size at each target complexity */ + /*!\brief ADAP_HMINS \n DESCRIPTION: Minimum cell size at each target complexity */ addPythonOption("ADAP_HMINS"); /*!\brief ADAP_HGRAD \n DESCRIPTION: Size gradation smoothing parameter */ addPythonOption("ADAP_HGRAD"); - /*!\brief ADAP_ITER \n DESCRIPTION: Hausdorff distance parameter for surface remeshing */ + /*!\brief ADAP_HAUSD \n DESCRIPTION: Hausdorff distance parameter for surface remeshing */ addPythonOption("ADAP_HAUSD"); - /*!\brief ADAP_ITER \n DESCRIPTION: A mesh adaptation option */ + /*!\brief ADAP_ANGLE \n DESCRIPTION: Sharp angle detection parameter for surface remeshing */ addPythonOption("ADAP_ANGLE"); /* END_CONFIG_OPTIONS */ From 79da503bd55b6da4e4d1cf3b8add58d643a6b052 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Tue, 4 Nov 2025 13:45:54 -0800 Subject: [PATCH 09/12] Fix call to SetPrimitive_Adapt() --- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index eec6ec19bc9a..f4343ba9d150 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -334,7 +334,7 @@ void CSinglezoneDriver::ComputeMetricField() { solver_flow->CompleteComms(geometry, config, MPI_QUANTITIES::SOLUTION); solver_flow->Preprocessing(geometry, solver, config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); - solver_flow->SetPrimitive_Adapt(geometry, config, solver_flow->GetNodes()); + solver_flow->SetPrimitive_Adapt(geometry, config); if (config->GetKind_Hessian_Method() == GREEN_GAUSS) { if(rank == MASTER_NODE) cout << "Computing Hessians using Green-Gauss." << endl; From b4c9765f1501e06ced62a6eafb8214e6686d7543 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Tue, 4 Nov 2025 14:17:27 -0800 Subject: [PATCH 10/12] Remove various unused variables and goal-oriented code paths --- Common/include/CConfig.hpp | 6 -- Common/src/CConfig.cpp | 12 ++- SU2_CFD/include/gradients/computeMetrics.hpp | 2 - SU2_CFD/include/solvers/CEulerSolver.hpp | 87 +++++++++----------- SU2_CFD/include/solvers/CIncEulerSolver.hpp | 73 ++++++++-------- SU2_CFD/include/solvers/CNEMOEulerSolver.hpp | 85 +++++++++---------- SU2_CFD/src/solvers/CSolver.cpp | 19 ++--- SU2_CFD/src/variables/CFlowVariable.cpp | 2 +- SU2_CFD/src/variables/CVariable.cpp | 2 +- 9 files changed, 127 insertions(+), 161 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index c2d7e854bd75..6f837fd0bc54 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -10166,12 +10166,6 @@ class CConfig { */ bool GetNormalize_Metric(void) const { return Normalize_Metric; } - /*! - * \brief Check if goal-oriented error estimation is being carried out - * \return TRUE<\code> if goal-oriented error estimation is taking place - */ - bool GetGoal_Oriented_Metric(void) const { return (nMetric_Sensor > 0) && (Metric_Sensor[0] == METRIC_SENSOR::GOAL); } - /*! * \brief Get the kind of method for computation of Hessians used for anisotropy. * \return Numerical method for computation of Hessians used for anisotropy. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 95a993c46eb6..956d6bcf4653 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5805,16 +5805,14 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Compute_Metric) { /*--- Check that config is valid for requested sensor ---*/ for (auto iSensor = 0; iSensor < nMetric_Sensor; iSensor++) { - /*--- If using GOAL, it must be the only sensor and the discrete adjoint must be used ---*/ + /*--- TODO: If using GOAL, it must be the only sensor and the discrete adjoint must be used ---*/ + /*--- For now, goal-oriented adaptation is unsupported ---*/ if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL) { - if (nMetric_Sensor != 1) - SU2_MPI::Error("Adaptation sensor GOAL cannot be used with other sensors.", CURRENT_FUNCTION); - if (!DiscreteAdjoint) - SU2_MPI::Error("Adaptation sensor GOAL can only be computed for MATH_PROBLEM = DISCRETE_ADJOINT.", CURRENT_FUNCTION); + SU2_MPI::Error("Adaptation sensor GOAL not yet supported.", CURRENT_FUNCTION); } if (Kind_Solver == MAIN_SOLVER::NEMO_EULER || Kind_Solver == MAIN_SOLVER::NEMO_NAVIER_STOKES) { - if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY || Metric_Sensor[iSensor] == METRIC_SENSOR::TOTAL_PRESSURE) + if (Metric_Sensor[iSensor] == METRIC_SENSOR::DENSITY || Metric_Sensor[iSensor] == METRIC_SENSOR::TOTAL_PRESSURE) SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for NEMO problems."), CURRENT_FUNCTION); } if (Kind_Solver != MAIN_SOLVER::NEMO_EULER && Kind_Solver != MAIN_SOLVER::NEMO_NAVIER_STOKES) { @@ -5822,7 +5820,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for non-NEMO problems."), CURRENT_FUNCTION); } if (Kind_Solver == MAIN_SOLVER::INC_EULER || Kind_Solver == MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver == MAIN_SOLVER::INC_RANS) { - if (Metric_Sensor[iSensor] == METRIC_SENSOR::GOAL || Metric_Sensor[iSensor] == METRIC_SENSOR::MACH) + if (Metric_Sensor[iSensor] == METRIC_SENSOR::MACH) SU2_MPI::Error(string("Adaptation sensor ") + GetMetric_SensorString(iSensor) + string(" not available for INC problems."), CURRENT_FUNCTION); } } diff --git a/SU2_CFD/include/gradients/computeMetrics.hpp b/SU2_CFD/include/gradients/computeMetrics.hpp index cb437fbd29ea..e1b90c17cc81 100644 --- a/SU2_CFD/include/gradients/computeMetrics.hpp +++ b/SU2_CFD/include/gradients/computeMetrics.hpp @@ -252,8 +252,6 @@ void normalizeMetrics(CGeometry& geometry, const CConfig& config, const unsigned long nPointDomain = geometry.GetnPointDomain(); - const bool goal = (config.GetGoal_Oriented_Metric()); - /*--- Constants defining normalization ---*/ const ScalarType p = config.GetMetric_Norm(); const ScalarType N = SU2_TYPE::GetValue(config.GetMetric_Complexity()); diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 0109d461015e..d8a4286d20cc 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -279,54 +279,49 @@ class CEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { - const auto nSensor = config->GetnMetric_Sensor(); - for (auto iSensor = 0; iSensor < nSensor; iSensor++) { - switch (config->GetMetric_Sensor(iSensor)) { - case METRIC_SENSOR::DENSITY: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetDensity(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::PRESSURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetPressure(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TOTAL_PRESSURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double Mach = sqrt(nodes->GetVelocity2(iPoint))/nodes->GetSoundSpeed(iPoint); - const su2double prim_var = nodes->GetPressure(iPoint)*pow((1+0.2*Mach*Mach), 3.5); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TEMPERATURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetTemperature(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::ENERGY: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetEnergy(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::MACH: - default: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - } + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::DENSITY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetDensity(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TOTAL_PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double Mach = sqrt(nodes->GetVelocity2(iPoint))/nodes->GetSoundSpeed(iPoint); + const su2double prim_var = nodes->GetPressure(iPoint)*pow((1+0.2*Mach*Mach), 3.5); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::MACH: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; } } - else { - // TODO: store variables which need Hessian computation for goal-oriented metric - } } /*! diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index c9fe7191d6ff..941f46c2eb52 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -92,47 +92,42 @@ class CIncEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { - const auto nSensor = config->GetnMetric_Sensor(); - for (auto iSensor = 0; iSensor < nSensor; iSensor++) { - switch (config->GetMetric_Sensor(iSensor)) { - case METRIC_SENSOR::DENSITY: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetDensity(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TOTAL_PRESSURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetPressure(iPoint) + 0.5*nodes->GetDensity(iPoint)*nodes->GetVelocity2(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TEMPERATURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetTemperature(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::ENERGY: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetEnergy(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::PRESSURE: - default: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetPressure(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - } + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::DENSITY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetDensity(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TOTAL_PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint) + 0.5*nodes->GetDensity(iPoint)*nodes->GetVelocity2(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::PRESSURE: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; } } - else { - // TODO: store variables which need Hessian computation for goal-oriented metric - } } /*! diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index 1f13399554c3..3b0dc0d8389a 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -201,53 +201,48 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseGetGoal_Oriented_Metric()) { - const auto nSensor = config->GetnMetric_Sensor(); - for (auto iSensor = 0; iSensor < nSensor; iSensor++) { - switch (config->GetMetric_Sensor(iSensor)) { - case METRIC_SENSOR::PRESSURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetPressure(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TEMPERATURE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetTemperature(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::TEMPERATURE_VE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetTemperature_ve(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::ENERGY: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = nodes->GetEnergy(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::ENERGY_VE: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = *nodes->GetEve(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - case METRIC_SENSOR::MACH: - default: - for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); - nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); - } - break; - } + const auto nSensor = config->GetnMetric_Sensor(); + for (auto iSensor = 0; iSensor < nSensor; iSensor++) { + switch (config->GetMetric_Sensor(iSensor)) { + case METRIC_SENSOR::PRESSURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetPressure(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::TEMPERATURE_VE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetTemperature_ve(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = nodes->GetEnergy(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::ENERGY_VE: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = *nodes->GetEve(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; + case METRIC_SENSOR::MACH: + default: + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + const su2double prim_var = sqrt(nodes->GetVelocity2(iPoint)) / nodes->GetSoundSpeed(iPoint); + nodes->SetPrimitive_Adapt(iPoint, iSensor, prim_var); + } + break; } } - else { - // TODO: store variables which need Hessian computation for goal-oriented metric - } } /*! diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 16bbfc1ea3b1..2e47e50fdaf0 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1381,11 +1381,11 @@ void CSolver::GetCommCountAndType(const CConfig* config, MPI_TYPE = COMM_TYPE_DOUBLE; break; case MPI_QUANTITIES::GRADIENT_ADAPT: - COUNT_PER_POINT = config->GetGoal_Oriented_Metric()? nVar*nDim : config->GetnMetric_Sensor()*nDim; + COUNT_PER_POINT = config->GetnMetric_Sensor()*nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; case MPI_QUANTITIES::HESSIAN: - COUNT_PER_POINT = config->GetGoal_Oriented_Metric()? nVar*nSymMat : config->GetnMetric_Sensor()*nSymMat; + COUNT_PER_POINT = config->GetnMetric_Sensor()*nSymMat; MPI_TYPE = COMM_TYPE_DOUBLE; break; default: @@ -2191,9 +2191,9 @@ void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, const CConfig *config } void CSolver::SetHessian_GG(CGeometry *geometry, const CConfig *config, short idxVel, const unsigned short Kind_Solver) { - const auto& solution = config->GetGoal_Oriented_Metric()? base_nodes->GetSolution() : base_nodes->GetPrimitive_Adapt(); + const auto& solution = base_nodes->GetPrimitive_Adapt(); auto& gradient = base_nodes->GetGradient_Adapt(); - auto nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + auto nHess = config->GetnMetric_Sensor(); computeGradientsGreenGauss(this, MPI_QUANTITIES::GRADIENT_ADAPT, PERIODIC_GRAD_ADAPT, *geometry, *config, solution, 0, nHess, idxVel, gradient); @@ -4391,14 +4391,10 @@ void CSolver::ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig /*--- TODO: - goal-oriented metric ---*/ /*--- - metric intersection ---*/ const unsigned long nPointDomain = geometry->GetnPointDomain(); + const unsigned short nSensor = config->GetnMetric_Sensor(); - const bool visc = (config->GetViscous()); - const bool turb = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); - - const bool goal = (config->GetGoal_Oriented_Metric()); const bool normalize = (config->GetNormalize_Metric()); - unsigned short nSensor = config->GetnMetric_Sensor(); const unsigned long time_iter = config->GetTimeIter(); const bool steady = (config->GetTime_Marching() == TIME_MARCHING::STEADY); @@ -4452,11 +4448,6 @@ void CSolver::AddMetrics(CSolver **solver, const CGeometry*geometry, const CConf const unsigned long nPointDomain = geometry->GetnPointDomain(); const unsigned short nSymMat = 3*(nDim-1); - const unsigned short nVarFlo = solver[FLOW_SOL]->GetnVar(); - const unsigned short nSensor = config->GetnMetric_Sensor(); - - const bool turb = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); - const bool goal = (config->GetGoal_Oriented_Metric()); const unsigned long time_iter = config->GetTimeIter(); const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || diff --git a/SU2_CFD/src/variables/CFlowVariable.cpp b/SU2_CFD/src/variables/CFlowVariable.cpp index d9350462ae19..2723543cee21 100644 --- a/SU2_CFD/src/variables/CFlowVariable.cpp +++ b/SU2_CFD/src/variables/CFlowVariable.cpp @@ -99,7 +99,7 @@ CFlowVariable::CFlowVariable(unsigned long npoint, unsigned long ndim, unsigned } if (config->GetCompute_Metric()) { - unsigned short nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + unsigned short nHess = config->GetnMetric_Sensor(); unsigned short nSymMat = 3 * (nDim - 1); Primitive_Adapt.resize(nPoint, nHess) = su2double(0.0); Metric.resize(nPoint, nSymMat) = 0.0; diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 8e7234b258ca..7c4ba2a70c16 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -83,7 +83,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva /*--- Gradient and Hessian for anisotropic metric ---*/ if (config->GetCompute_Metric()) { - unsigned short nHess = config->GetGoal_Oriented_Metric()? nVar : config->GetnMetric_Sensor(); + unsigned short nHess = config->GetnMetric_Sensor(); Gradient_Adapt.resize(nPoint,nHess,nDim,0.0); Hessian.resize(nPoint,nHess,nSymMat,0.0); } From 76f7a9d15a911c810894e032f96b5b04fdc88a17 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Tue, 4 Nov 2025 14:41:43 -0800 Subject: [PATCH 11/12] Address CodeQL complaints --- SU2_CFD/include/gradients/computeMetrics.hpp | 5 ++--- SU2_CFD/src/output/CFlowOutput.cpp | 1 - SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp | 1 - 3 files changed, 2 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/gradients/computeMetrics.hpp b/SU2_CFD/include/gradients/computeMetrics.hpp index e1b90c17cc81..8917d580973b 100644 --- a/SU2_CFD/include/gradients/computeMetrics.hpp +++ b/SU2_CFD/include/gradients/computeMetrics.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/linear_algebra/blas_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" @@ -171,7 +172,7 @@ void setPositiveDefiniteMetrics(CGeometry& geometry, const CConfig& config, /*--- Make positive definite by taking absolute value of eigenvalues ---*/ /*--- Handle NaN and very small values that could cause numerical issues ---*/ for (auto iDim = 0; iDim < nDim; iDim++) { - if (EigVal[iDim] != EigVal[iDim]) { + if (std::isnan(EigVal[iDim])) { /*--- NaN detected, set to small positive value ---*/ EigVal[iDim] = eps; } else { @@ -268,8 +269,6 @@ void normalizeMetrics(CGeometry& geometry, const CConfig& config, ScalarType A[nDim][nDim], EigVec[nDim][nDim], EigVal[nDim], work[nDim]; for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { - auto nodes = geometry.nodes; - /*--- Decompose metric ---*/ Tensor::get(metric, iPoint, iSensor, A, nDim); CBlasStructure::EigenDecomposition(A, EigVec, EigVal, nDim, work); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index be5647734a0b..f6c7c38735d4 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -4132,7 +4132,6 @@ void CFlowOutput::LoadMeshAdaptationOutputs(const CConfig* config, const CSolver unsigned long iPoint) { const auto* Node_Flow = solver[FLOW_SOL]->GetNodes(); - const auto* Node_Geo = geometry->nodes; if(config->GetCompute_Metric()) { // Common metric components for both 2D and 3D SetVolumeOutputValue("METRIC_XX", iPoint, Node_Flow->GetMetric(iPoint, 0)); diff --git a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp index 14fd3e1a3fa2..da99b870fe4b 100644 --- a/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp +++ b/SU2_CFD/src/output/filewriter/CParaviewXMLFileWriter.cpp @@ -384,7 +384,6 @@ void CParaviewXMLFileWriter::WriteData(string val_filename){ /*--- Resize buffer to accommodate tensor data ---*/ dataBufferFloat.resize(myPoint*NTENSOR); - float val = 0.0; for (iPoint = 0; iPoint < myPoint; iPoint++) { dataBufferFloat[iPoint*NTENSOR + 0] = (float)dataSorter->GetData(VarCounter+0,iPoint); // XX dataBufferFloat[iPoint*NTENSOR + 1] = (float)dataSorter->GetData(VarCounter+2,iPoint); // YY From 6e4c7606ef4f7ac321630f9a86b83641afb08da6 Mon Sep 17 00:00:00 2001 From: bmunguia Date: Tue, 4 Nov 2025 18:15:45 -0800 Subject: [PATCH 12/12] Bug fix --- SU2_CFD/src/solvers/CSolver.cpp | 47 +++++++++++++++------------------ 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 2e47e50fdaf0..58fc19bd293f 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4407,33 +4407,28 @@ void CSolver::ComputeMetric(CSolver **solver, CGeometry *geometry, const CConfig vector integrals; for (auto iSensor = 0u; iSensor < nSensor; ++iSensor) { SU2_OMP_MASTER - - if (!goal) { - /*--- Make the Hessian eigenvalues positive definite, and add to the metric tensor ---*/ - auto& hessians = base_nodes->GetHessian(); - setPositiveDefiniteMetrics(*geometry, *config, iSensor, hessians); - AddMetrics(solver, geometry, config, iSensor); - - /*--- Integrate metric field on the last iteration (the end of the simulation if steady) ---*/ - auto& metrics = base_nodes->GetMetric(); - double integral = 0.0; - if (is_last_iter) - integral = integrateMetrics(*geometry, *config, iSensor, metrics); - - /*--- Normalize the metric field for steady simulations, or if requested for unsteady ---*/ - if (steady || (normalize && is_last_iter)) - normalizeMetrics(*geometry, *config, iSensor, integral, metrics); - - /*--- Store the integral to be written ---*/ - if (is_last_iter) { - integrals.push_back(integral); - if (rank == MASTER_NODE) { - cout << "Global metric normalization integral for sensor "; - cout << config->GetMetric_SensorString(iSensor) << ": " << integral << endl; - } + /*--- Make the Hessian eigenvalues positive definite, and add to the metric tensor ---*/ + auto& hessians = base_nodes->GetHessian(); + setPositiveDefiniteMetrics(*geometry, *config, iSensor, hessians); + AddMetrics(solver, geometry, config, iSensor); + + /*--- Integrate metric field on the last iteration (the end of the simulation if steady) ---*/ + auto& metrics = base_nodes->GetMetric(); + double integral = 0.0; + if (is_last_iter) + integral = integrateMetrics(*geometry, *config, iSensor, metrics); + + /*--- Normalize the metric field for steady simulations, or if requested for unsteady ---*/ + if (steady || (normalize && is_last_iter)) + normalizeMetrics(*geometry, *config, iSensor, integral, metrics); + + /*--- Store the integral to be written ---*/ + if (is_last_iter) { + integrals.push_back(integral); + if (rank == MASTER_NODE) { + cout << "Global metric normalization integral for sensor "; + cout << config->GetMetric_SensorString(iSensor) << ": " << integral << endl; } - } else { - SU2_MPI::Error("Goal-oriented metric not currently implemented.", CURRENT_FUNCTION); } END_SU2_OMP_MASTER SU2_OMP_BARRIER