ENH: changes in advectionDiffusionPatchDistMethod

- advectionDiffusion is frequently used within optimisation loops since
  it is differentiable. In shape optimisation, the re-computation of
  mesh distances is performed at the very beginning of a new
  optimisation cycle, due to inheriting from MeshObject. If the mesh
  quality is poor enough, the advectionDiffusion PDE might diverge and
  crash the run, before the problematic mesh is written to files for
  inspection. The default behaviour now is to check the mesh before
  solving the advectionDiffusion PDE and write the mesh points if some
  mesh check fails.
- fvOptions can now be included in advectionDiffusion (necessary for
  topology optimisation of turbulent flows for models that include the
  distance field)
- Minor changes in the numerical treatment of the diffusion term, to
  enhance stability

src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C

@@ -7,6 +7,8 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2017 OpenFOAM Foundation
     Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020-2023 PCOpt/NTUA
+    Copyright (C) 2020-2023 FOSS GP
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +35,7 @@ License
 #include "fvmDiv.H"
 #include "fvmLaplacian.H"
 #include "fvmSup.H"
+#include "fvOptions.H"
 #include "addToRunTimeSelectionTable.H"
 
 #include "fixedValueFvPatchFields.H"
@@ -75,7 +78,9 @@ Foam::patchDistMethods::advectionDiffusion::advectionDiffusion
     epsilon_(coeffs_.getOrDefault<scalar>("epsilon", 0.1)),
     tolerance_(coeffs_.getOrDefault<scalar>("tolerance", 1e-3)),
     maxIter_(coeffs_.getOrDefault<int>("maxIter", 10)),
-    predicted_(false)
+    predicted_(false),
+    checkAndWriteMesh_
+        (coeffs_.lookupOrDefault<bool>("checkAndWriteMesh", true))
 {}
 
 
@@ -99,6 +104,33 @@ bool Foam::patchDistMethods::advectionDiffusion::correct
         predicted_ = true;
     }
 
+    // If the mesh has become invalid, trying to solve the eikonal equation
+    // might lead to crashing and we wont get a chance to see the failed mesh.
+    // Write the mesh points here
+    if (checkAndWriteMesh_)
+    {
+        DebugInfo
+            << "Checking mesh in advectionDiffusion " << endl;
+        if (mesh_.checkMesh(true))
+        {
+            pointIOField points
+            (
+                IOobject
+                (
+                   "points",
+                    mesh_.pointsInstance(),
+                    mesh_.meshSubDir,
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh_.points()
+            );
+            points.write();
+        }
+    }
+
     volVectorField ny
     (
         IOobject
@@ -123,6 +155,23 @@ bool Foam::patchDistMethods::advectionDiffusion::correct
         nybf[patchi] == -patches[patchi].nf();
     }
 
+    // Scaling dimensions, to be used with fvOptions
+    volScalarField scaleDims
+    (
+        IOobject
+        (
+            "scaleDims",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        mesh_,
+        dimensionedScalar("scaleDims", dimTime/dimLength, scalar(1))
+    );
+
+    fv::options& fvOptions(fv::options::New(this->mesh_));
     int iter = 0;
     scalar initialResidual = 0;
 
@@ -139,14 +188,17 @@ bool Foam::patchDistMethods::advectionDiffusion::correct
         fvScalarMatrix yEqn
         (
             fvm::div(yPhi, y)
-          - fvm::Sp(fvc::div(yPhi), y)
+          + fvm::SuSp(-fvc::div(yPhi), y)
           - epsilon_*y*fvm::laplacian(y)
          ==
             dimensionedScalar("1", dimless, 1.0)
+          + fvOptions(scaleDims, y)
         );
 
         yEqn.relax();
+        fvOptions.constrain(yEqn);
         initialResidual = yEqn.solve().initialResidual();
+        fvOptions.correct(y);
 
     } while (initialResidual > tolerance_ && ++iter < maxIter_);
 

src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.H

@@ -6,6 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2016 OpenFOAM Foundation
+    Copyright (C) 2020-2023 PCOpt/NTUA
+    Copyright (C) 2020-2023 FOSS GP
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -206,6 +208,18 @@ class advectionDiffusion
         //- Flag to indicate the predictor step has been executed
         bool predicted_;
 
+        //- Check mesh before computing the distance field and write the mesh
+        //- points if at least one check fails.
+        //  advectionDiffusion is frequently used within optimisation loops
+        //  since it is differentiable. In shape optimisation, the
+        //  re-computation of mesh distances is performed at the very beginning
+        //  of a new optimisation cycle, due to inheriting from MeshObject. If
+        //  the mesh quality is poor enough, the advectionDiffusion PDE might
+        //  diverge and crash the run, before the problematic mesh is written
+        //  to files for inspection. By setting this flag to true, the mesh is
+        //  written to files before the solution of the advectionDiffusion PDE,
+        //  if some mesh check fails.
+        bool checkAndWriteMesh_;
 
     // Private Member Functions