From bf18fb758c629b208b1ab811162a24f58910ec22 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Thu, 23 Nov 2023 14:48:11 +0200
Subject: [PATCH] ENH: nullSpace will now perturb the design variables

by a small amount, if all of them lay on the lower or upper bounds at
the beginning of the optimisation, to avoid singular matrices when
computing the update of the design variables.
---
 .../updateMethod/nullSpace/nullSpace.C        | 59 +++++++++++++++++++
 .../updateMethod/nullSpace/nullSpace.H        |  4 ++
 2 files changed, 63 insertions(+)

diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.C
index 7377f90078..82a4207ab0 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.C
@@ -971,6 +971,7 @@ Foam::nullSpace::nullSpace
         coeffsDict(type).
             getOrDefault<scalar>("violatedConstraintsThreshold", 1.e-3)
     ),
+    epsPerturb_(coeffsDict(type).  getOrDefault<scalar>("perturbation", 1.e-2)),
     aJ_(coeffsDict(type).getOrDefault<scalar>("aJ", 1)),
     aC_(coeffsDict(type).getOrDefault<scalar>("aC", 1)),
     adaptiveStep_(coeffsDict(type).getOrDefault<bool>("adaptiveStep", true)),
@@ -1003,6 +1004,64 @@ Foam::nullSpace::nullSpace
             << exit(FatalError);
     }
 
+    // Sanity checks for the bounds of the design variables.
+    // If all design variables start from either their lower or upper bound
+    // and there is at least one active flow-related constraint, the systems
+    // providing the update of the desgin varaibles become singular.
+    // Issue a warning and change the initial values of the design variables
+    // by a small ammount, to kick start the optimisation
+    if (designVars_->lowerBounds() && designVars_->upperBounds())
+    {
+        const scalarField& lowerBounds = designVars_->lowerBoundsRef();
+        const scalarField& upperBounds = designVars_->upperBoundsRef();
+        scalarField& designVars = designVars_();
+        boolList isOnLowerBound(designVars.size(), false);
+        boolList isOnUpperBound(designVars.size(), false);
+        bool existsNonBoundVar = false;
+        for (const label activeVarI : designVars_->activeDesignVariables())
+        {
+            isOnLowerBound[activeVarI] =
+                mag(designVars[activeVarI] - lowerBounds[activeVarI]) < SMALL;
+            isOnUpperBound[activeVarI] =
+                mag(designVars[activeVarI] - upperBounds[activeVarI]) < SMALL;
+            existsNonBoundVar =
+                existsNonBoundVar
+            || (!isOnLowerBound[activeVarI] && !isOnUpperBound[activeVarI]);
+        }
+        if (globalSum_)
+        {
+            reduce(existsNonBoundVar, orOp<bool>());
+        }
+        if (!existsNonBoundVar)
+        {
+            WarningInFunction
+                << "All design variables lay on their bound values." << nl
+                << tab << "This will lead to singular matrices in nullSpace"
+                << nl
+                << tab << "Perturbing the design variables by "
+                << epsPerturb_ << "*(upperBound - lowerBounds)"
+                << nl
+                << endl;
+            scalarField update(designVars.size(), Zero);
+            for (const label activeVarI : designVars_->activeDesignVariables())
+            {
+                if (isOnLowerBound[activeVarI])
+                {
+                    update[activeVarI] =
+                        epsPerturb_
+                       *(upperBounds[activeVarI] - lowerBounds[activeVarI]);
+                }
+                if (isOnUpperBound[activeVarI])
+                {
+                    update[activeVarI] =
+                      - epsPerturb_
+                       *(upperBounds[activeVarI] - lowerBounds[activeVarI]);
+                }
+            }
+            designVars_->update(update);
+        }
+    }
+
     // Read-in aJ in restarts
     dictionary::readIfPresent("aJ", aJ_);
 }
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.H
index 951adf3f42..294f66fb1c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.H
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/updateMethod/nullSpace/nullSpace.H
@@ -159,6 +159,10 @@ protected:
             //  constraints
             scalar epsConstr_;
 
+            //- Value to perturb the design variables by, if all of them
+            //- lay on their bounds at the beginning of the optimisation
+            scalar epsPerturb_;
+
             //- Multiplier of the null space update
             scalar aJ_;