From 59bf69b92e7cb16a56f031c9502a764fa26b539e Mon Sep 17 00:00:00 2001 From: Vaggelis Papoutsis Date: Wed, 19 Jul 2023 12:10:14 +0300 Subject: [PATCH] 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 --- .../advectionDiffusionPatchDistMethod.C | 56 ++++++++++++++++++- .../advectionDiffusionPatchDistMethod.H | 14 +++++ 2 files changed, 68 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C index 283cbb4331..915b62d003 100644 --- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C +++ b/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("epsilon", 0.1)), tolerance_(coeffs_.getOrDefault("tolerance", 1e-3)), maxIter_(coeffs_.getOrDefault("maxIter", 10)), - predicted_(false) + predicted_(false), + checkAndWriteMesh_ + (coeffs_.lookupOrDefault("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_); diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.H b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.H index 6e5d4d443e..c8201247bf 100644 --- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.H +++ b/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