Anisotropic mesh adaptation part 1: feature-based metric tensor

Common/include/CConfig.hpp

@@ -1295,6 +1295,20 @@ 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. */
+  string* Metric_Sensor;                   /*!< \brief Sensors to use for adaptation (first entry is normalized, rest are Hessian-only). */
+
+  unsigned short Metric_Norm;              /*!< \brief Lp-norm for mesh adaptation */
+  unsigned long Metric_Complexity;         /*!< \brief Constraint mesh complexity */
+  unsigned short nAdapt_Time_Subinterval;  /*!< \brief Number of unsteady time sub-intervals for adaptation. */
+  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.
@@ -10244,4 +10258,91 @@ class CConfig {
    */
   const FluidFlamelet_ParsedOptions& GetFlameletParsedOptions() const { return flamelet_ParsedOptions; }
 
+  /*!
+   * \brief Check if error estimation is being carried out
+   * \return <code>TRUE<\code> if error estimation is taking place
+  */
+  bool GetCompute_Metric(void) const { return Compute_Metric; }
+
+  /*!
+   * \brief Check if metric tensor normalization is being carried out
+   * \return <code>TRUE<\code> if metric normalization is taking place
+  */
+  bool GetNormalize_Metric(void) const { return Normalize_Metric; }
+
+  /*!
+   * \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 complete array of metric sensor names
+   * \return Array of sensor names
+   */
+  string* GetMetric_Sensor() const {
+    return Metric_Sensor;
+  }
+
+  /*!
+   * \brief Get metric sensor name by index
+   * \param[in] iSens - Index of the sensor
+   * \return Sensor name string
+   */
+  string GetMetric_Sensor(unsigned short iSens) const {
+    if (iSens >= nMetric_Sensor)
+      SU2_MPI::Error("Sensor index out of range.", CURRENT_FUNCTION);
+    return Metric_Sensor[iSens];
+  }
+
+  /*!
+   * \brief Get the complete list of metric sensor names
+   * \return Vector of sensor name strings
+   */
+  vector<string> GetMetric_SensorList() const {
+    return vector<string>(Metric_Sensor, Metric_Sensor + nMetric_Sensor);
+  }
+
+  /*!
+   * \brief Get number of adaptation sensors
+   * \return Number of sensors
+   */
+  unsigned short GetnMetric_Sensor() 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; }
+
+  /*!
+   * \brief Get number of unsteady adaptation sub-intervals
+   * \note Currently only one sub-interval supported
+   * \return Number of unsteady adaptation sub-intervals
+   */
+  unsigned long GetnAdapt_Time_Subinterval(void) const { return nAdapt_Time_Subinterval; }
+
 };

Common/include/option_structure.hpp

@@ -2686,6 +2686,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). */
 };
 
 /*!
@@ -2717,6 +2719,9 @@ enum class MPI_QUANTITIES {
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
+  SENSOR_ADAPT         ,  /*!< \brief Sensors for anisotropic sizing metric tensor. */
+  GRADIENT_ADAPT       ,  /*!< \brief Gradient vectors for anisotropic sizing metric tensor. */
+  HESSIAN              ,  /*!< \brief Hessian tensors for anisotropic sizing metric tensor. */
 };
 
 /*!
@@ -2874,6 +2879,15 @@ static const MapType<std::string, DEFORM_KIND> Deform_Kind_Map = {
   MakePair("RBF",       DEFORM_KIND::RBF)
 };
 
+/*!
+ * \brief Type of sensor for anisotropic metrics.
+ */
+enum class SensorType : unsigned char {
+  PRIMITIVE, /*!< \brief Value read directly from the primitive variable array. */
+  DERIVED,   /*!< \brief Officially-supported computed quantity (e.g. Mach number). */
+  CUSTOM,    /*!< \brief User-defined sensor populated externally via the Python wrapper. */
+};
+
 
 #undef MakePair
 /* END_CONFIG_ENUMS */

Common/src/CConfig.cpp

@@ -3083,6 +3083,60 @@
   /*!\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 an error estimate */
+  addBoolOption("COMPUTE_METRIC", Compute_Metric, false);
+  /*!\brief NORMALIZE_METRIC \n DESCRIPTION: Normalize the metric tensor */
+  addBoolOption("NORMALIZE_METRIC", Normalize_Metric, false);
+  /*!\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 metric field */
+  addStringListOption("METRIC_SENSOR", nMetric_Sensor, Metric_Sensor);
+  /*!\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 ADAP_TIME_SUBINTERVAL \n DESCRIPTION: Number of time subintervals in unsteady mesh adaptation */
+  addUnsignedShortOption("ADAP_TIME_SUBINTERVAL", nAdapt_Time_Subinterval, 1);
+
+  /*!\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 METRIC_HGRAD \n DESCRIPTION: Size gradation smoothing parameter */
+  addPythonOption("METRIC_HGRAD");
+
+  /*!\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_HMAXS \n DESCRIPTION: Maximum cell size at each target complexity */
+  addPythonOption("ADAP_HMAXS");
+  /*!\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_HAUSD \n DESCRIPTION: Hausdorff distance parameter for surface remeshing */
+  addPythonOption("ADAP_HAUSD");
+  /*!\brief ADAP_ANGLE \n DESCRIPTION: Sharp angle detection parameter for surface remeshing */
+  addPythonOption("ADAP_ANGLE");
+
   /* END_CONFIG_OPTIONS */
 
 }
@@ -5713,6 +5767,30 @@
     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 (unsigned short iSensor = 0; iSensor < nMetric_Sensor; iSensor++) {
+      const string& sensor_name = Metric_Sensor[iSensor];
+      /*--- If using GOAL, it must be the only sensor and the discrete adjoint must be used ---*/
+      /*--- TODO: goal-oriented adaptation ---*/
+      if (sensor_name == "GOAL") {
+        SU2_MPI::Error("Adaptation sensor GOAL not yet supported.", 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);
+    }
+
+    /*--- Make sure only using single adaptation sub-interval for steady problems ---*/
+    if(TimeMarching == TIME_MARCHING::STEADY)
+      nAdapt_Time_Subinterval = 1;
+    if (nAdapt_Time_Subinterval != 1)
+      SU2_MPI::Error("Adaptation sub-intervals not yet supported. Set ADAP_TIME_SUBINTERVAL = 1 or remove from config.", CURRENT_FUNCTION);
+  }
+
 }
 
 void CConfig::SetMarkers(SU2_COMPONENT val_software) {
@@ -7914,6 +7992,34 @@
       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 (unsigned short iSensor = 0; iSensor < nMetric_Sensor; iSensor++) {
+        cout << Metric_Sensor[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;
+        if (TimeMarching != TIME_MARCHING::STEADY) {
+          cout << "  Unsteady adaptation sub-intervals: " << nAdapt_Time_Subinterval << endl;
+          cout << "  Target space-time complexity: " << Metric_Complexity * nAdapt_Time_Subinterval << 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,

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -515,6 +515,23 @@ class CDriverBase {
    */
   map<string, unsigned short> GetPrimitiveIndices() const;
 
+  /*!
+   * \brief Get the local index of a named metric sensor in the flow solver.
+   * Used by Python custom sensor registries to cache indices before the run loop.
+   * \param[in] sensor_name - Name as listed in METRIC_SENSOR config.
+   * \return Sensor index, or -1 if not found.
+   */
+  short GetMetricSensorIndex(const std::string& sensor_name) const;
+
+  /*!
+   * \brief Get a read/write view of adaptation sensor values on all mesh nodes of the flow solver.
+   * \warning Adaptation sensors are only available for flow solvers with metric sensors configured.
+   */
+  inline CPyWrapperMatrixView AdaptSensors() {
+    auto* solver = GetSolverAndCheckMarker(FLOW_SOL);
+    return CPyWrapperMatrixView(const_cast<su2activematrix&>(solver->GetNodes()->GetSensor_Adapt()), "AdaptSensors", false);
+  }
+
   /*!
    * \brief Get a read/write view of the current primitive variables on all mesh nodes of the flow solver.
    * \warning Primitive variables are only available for flow solvers.

SU2_CFD/include/drivers/CSinglezoneDriver.hpp

@@ -110,4 +110,10 @@ class CSinglezoneDriver : public CDriver {
    */
   bool Monitor(unsigned long TimeIter) override;
 
+  /*!
+   * \brief Perform all steps to compute the metric tensor.
+   * \param[in] restartMetric - Whether this is the initial sub-interval metric computation for an unsteady restart.
+   */
+  virtual void ComputeMetricField(bool restartMetric = false);
+
 };