Anisotropic mesh adaptation part 1: feature-based metric tensor

[bmunguia]

Proposed Changes

This PR is part of a joint effort with Alberto Perlini (@xla27) and Andrea Rausa (@rois1995) to implement anisotropic mesh adaptation in SU2. In our SU2 conference talks this year, I presented updates to the metric tensor computation in my fork for unsteady adaptation, and Alberto presented updates to the Python adaptation loop to use the open-source MMG remesher, as well as results for incompressible Euler/NS and NEMO.

Part 1 will focus on my contributions in SU2's post-processing/output to calculate metric tensor fields for feature-based adaptation. Part 2 will focus on the implementation of the adaptation loop in Python, including conversion between SU2 restarts/meshes and the .meshb format.

For the metric tensor field computation, I have implemented:

• CVariable containers to store adaptation sensor variable/gradient/Hessian/metric tensor
  • We only store the lower-triangular elements of the Hessian/metric, so the size of the array at each point is $nSymMat = 3*(nDim-1)$
• CSolver function SetPrimitive_Adapt() to store the chosen adaptation sensor
• Green-Gauss Hessian calculation in computeGradientsGreenGauss.hpp which calculates the Hessian field for the chosen sensor and stores the lower triangle
• Integration of the Hessian (symmetric positive definite by construction) determinant and Lp-norm normalization (see references) in computeMetrics.hpp
• Integration of the Hessian field in time for unsteady adaptation
• Output of metric tensor for different flow solvers, and updates to Paraview XML filewriter to allow for proper parsing of symmetric tensors

The following relevant config options have been added:

• COMPUTE_METRIC (bool) - whether to calculate the sizing metric in post-processing (default False)
• NORMALIZE_METRIC (bool) - whether to perform the Lp-norm normalization of Hessians (default True)
• NUM_METHOD_HESS (str) - gradient method for Hessian calculation (only GG for now)
• METRIC_SENSOR (str) - field to use for sizing metric
• METRIC_NORM (unsigned short) - choice of norm for Lp-norm normalization
• METRIC_COMPLEXITY (unsigned long) - target mesh size
• METRIC_HMAX (float) - maximum cell size (constraint on minimum metric tensor eigenvalue)
• METRIC_HMIN (float) - minimum cell size (constraint on maximum metric tensor eigenvalue)
• METRIC_ARMAX (float) - maximum cell aspect ratio (constraint on ratio of hmax to hmin)

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

References

1. Loseille and Alauzet. “Continuous mesh framework part I: well-posed continuous interpolation error,” 2011.
2. Alauzet and Loseille. “A decade of progress on anisotropic mesh adaptation for computational fluid dynamics,” 2016.
3. Perlini. "Mesh adaptation in SU2: current status and perspectives," 2025.
4. Munguía. "Towards a deep learning framework for time-accurate anisotropic mesh adaptation," 2025.

[pcarruscag]

Is the mesh adaptation code itself open source?

[bmunguia]

If you mean MMG, it's on Github https://github.com/MmgTools/mmg

If you mean the SU2-MMG interface, it's in Alberto's fork https://github.com/xla27/SU2/tree/feat_adap_mmg

[pcarruscag]

Excellent. Next question, what do you need from the python wrapper to do all of this outside of the SU2 core?

[bmunguia]

The python script basically just coordinates calls to SU2 and MMG in 2 nested loops, the inner loop being a series of adaptations at a given target mesh size and the outer loop being over a list of sizes. It also handles file conversions since MMG uses libmeshb for both the mesh and metric field.

For the file conversion, Alberto has a straight python implementation that I think only handles tris/tets, and I have a separate repo with a C++ implementation/python API that I used for my results that I presented in Varenna https://github.com/bmunguia/limon

[pcarruscag]

Cool, the functionality would be more flexible if the sensor was also configurable via python, what you have here has a list of hardcoded possible variables, in python it would be very easy to target other solvers, create custom mixes of variables etc.
So maybe a helper to compute a generic gradient/hessian in c++?

[bigfooted]

Hi Brian, what's the status?

[bmunguia]

Hey Nijso, sorry for the delay. I've been a bit swamped lately with a deadline and the holidays so this has been on the backburner, but I'll look at it again this week or next. I know a lot of people have been asking about adaptation on slack so I want to get this in develop ASAP.

I still need to implement Pedro's suggestions to use the python wrapper to make the choice of metric sensor more flexible, then I just need to add an example (or maybe a few for different solvers) to test the implementation.

[bigfooted]

Hi @bmunguia, any update on this?

[bmunguia]

Sorry for the massive delay. I made the implementation a bit more flexible by resolving the indices using CPrimitiveIndices rather than hardcoding per-solver. So now if you specify a list of sensors

COMPUTE_METRIC= YES
NORMALIZE_METRIC= YES
METRIC_SENSOR= DENSITY, PRESSURE, TEMPERATURE
METRIC_NORM= 2
WRT_RESTART_COMPACT= NO
VOLUME_OUTPUT= (COORDINATES, SOLUTION, PRIMITIVE, SENSOR_GRADIENT, SENSOR_HESSIAN, METRIC)

it should set them all as sensors

void CSolver::SetPrimitive_Adapt(CGeometry *geometry, const CConfig *config) {
  const auto nSensors = MetricSensorIndices.size();

  /*--- Copy each resolved sensor variable into Sensor_Adapt ---*/
  for (size_t iSensor = 0; iSensor < nSensors; iSensor++) {
    const auto var_idx = MetricSensorIndices[iSensor];

    SU2_OMP_FOR_STAT(omp_chunk_size)
    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
      const su2double prim_var = base_nodes->GetPrimitive(iPoint, var_idx);
      base_nodes->SetSensor_Adapt(iPoint, iSensor, prim_var);
    }
    END_SU2_OMP_FOR
  }
}

and calculate gradients/Hessians for all of them, as well as a Lp-norm normalized metric field for the first entry. I still need to add the python wrapper bit so we can compute metrics of derived quantities like Mach number. This is currently my top priority since I actually need the Mach metric soon for my own work. I'm hoping to get it implemented today or tomorrow.

[bmunguia]

The overhead of spawning the python wrapper just for a Mach Hessian is too much for my use case, so I implemented a compromise of the Polimi implementation and Pedro's request for a more flexible python wrapper-based implementation using 3 sensor types:

• Primitive: directly extracted from CVariable via PrimitiveNameToIndexMap() (e.g. DENSITY, PRESSURE)
• Derived: a limited number of supported derived quantities (e.g. MACH, VELOCITY)
• Custom: computed via the python wrapper (e.g. TOTAL_PRESSURE in my example)

[bigfooted]

Hi @bmunguia any updates on this? Can you resolve the conflicts, then we can finalize this implementation.

[bmunguia]

Hi @bigfooted it was more or less ready ~6 months ago when I addressed all the comments but I guess it slipped through the cracks. Anyway, I merged yesterday and made some small changes to how custom adaptation sensors are handled. Let me know if there's anything you want cleaned up or handled differently.

Once this is finalized, I can work with the Polimi students on updating their SU2-MMG interface.