multiphaseEuler: Prevent sigFpe errors from blending method NaN handling

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/blendingMethod/blendingMethod.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-Foam::scalar Foam::blendingMethod::readParameter
+Foam::blendingParameter Foam::blendingMethod::readParameter
 (
     const word& name,
     const dictionary& dict,
@@ -51,28 +51,25 @@ Foam::scalar Foam::blendingMethod::readParameter
 
         if (allowNone && t.isWord() && t.wordToken() == "none")
         {
-            return NaN;
+            return {false, NaN};
         }
 
         if (t.isNumber())
         {
             forAll(bounds, i)
             {
-                if (!std::isnan(bounds[i]))
+                const label s = i == 0 ? -1 : +1;
-                {
-                    const label s = i == 0 ? -1 : +1;
 
-                    if (s*t.number() > s*bounds[i])
+                if (s*t.number() > s*bounds[i])
-                    {
+                {
-                        FatalErrorInFunction
+                    FatalErrorInFunction
-                            << "Blending parameter " << name << " is "
+                        << "Blending parameter " << name << " is "
-                            << (i == 0 ? "less" : "greater") << " than "
+                        << (i == 0 ? "less" : "greater") << " than "
-                            << bounds[i] << exit(FatalError);
+                        << bounds[i] << exit(FatalError);
-                    }
                 }
             }
 
-            return t.number();
+            return {true, t.number()};
         }
 
         FatalIOErrorInFunction(dict)
@@ -80,11 +77,11 @@ Foam::scalar Foam::blendingMethod::readParameter
             << t.info() << exit(FatalIOError);
     }
 
-    return NaN;
+    return {false, NaN};
 }
 
 
-Foam::Pair<Foam::scalar> Foam::blendingMethod::readParameters
+Foam::Pair<Foam::blendingParameter> Foam::blendingMethod::readParameters
 (
     const word& name,
     const dictionary& dict,
@@ -96,7 +93,7 @@ Foam::Pair<Foam::scalar> Foam::blendingMethod::readParameters
     const word name1 = IOobject::groupName(name, interface.phase1().name());
     const word name2 = IOobject::groupName(name, interface.phase2().name());
     return
-        Pair<scalar>
+        Pair<blendingParameter>
         (
             readParameter(name1, dict, bounds, allowNone),
             readParameter(name2, dict, bounds, allowNone)
@@ -104,12 +101,6 @@ Foam::Pair<Foam::scalar> Foam::blendingMethod::readParameters
 }
 
 
-bool Foam::blendingMethod::isParameter(const scalar parameter)
-{
-    return !std::isnan(parameter);
-}
-
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::blendingMethod::constant
@@ -156,7 +147,7 @@ Foam::tmp<Foam::volScalarField> Foam::blendingMethod::parameter
 (
     const UPtrList<const volScalarField>& alphas,
     const label set,
-    const Pair<scalar>& parameters
+    const Pair<blendingParameter>& parameters
 ) const
 {
     tmp<volScalarField> talphaParameter = constant(alphas, 0);
@@ -167,7 +158,7 @@ Foam::tmp<Foam::volScalarField> Foam::blendingMethod::parameter
         {
             talphaParameter.ref() +=
                 max(iter().residualAlpha(), alphas[iter().index()])
-               *parameters[iter.index()];
+               *parameters[iter.index()].value;
         }
     }
 

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/blendingMethod/blendingMethod.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,20 @@ SourceFiles
 namespace Foam
 {
 
+/*---------------------------------------------------------------------------*\
+                     Struct blendingParameter Declaration
+\*---------------------------------------------------------------------------*/
+
+struct blendingParameter
+{
+    //- Is this parameter valid?
+    bool valid;
+
+    //- The parameter value. Could be templated if needed.
+    scalar value;
+};
+
+
 /*---------------------------------------------------------------------------*\
                        Class blendingMethod Declaration
 \*---------------------------------------------------------------------------*/
@@ -62,7 +76,7 @@ protected:
     // Protected Static Member Functions
 
         //- Read a parameter and check it lies within specified bounds
-        static scalar readParameter
+        static blendingParameter readParameter
         (
             const word& name,
             const dictionary& dict,
@@ -71,7 +85,7 @@ protected:
         );
 
         //- Read a parameter for each phase in the interface
-        static Pair<scalar> readParameters
+        static Pair<blendingParameter> readParameters
         (
             const word& name,
             const dictionary& dict,
@@ -80,9 +94,6 @@ protected:
             const bool allowNone
         );
 
-        //- Test if the given scalar is a valid parameter
-        static bool isParameter(const scalar parameter);
-
 
     // Protected Member Functions
 
@@ -106,7 +117,7 @@ protected:
         (
             const UPtrList<const volScalarField>& alphas,
             const label set,
-            const Pair<scalar>& parameters
+            const Pair<blendingParameter>& parameters
         ) const;
 
         //- Return the coordinate of the blending function

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ Foam::tmp<Foam::volScalarField> Foam::blendingMethods::hyperbolic::fContinuous
 {
     tmp<volScalarField> x = this->x(alphas, phaseSet, systemSet);
     tmp<volScalarField> a = parameter(alphas, phaseSet, minContinuousAlpha_);
-    return (1 + tanh((4/transitionAlphaScale_)*(x - a)))/2;
+    return (1 + tanh((4/transitionAlphaScale_.value)*(x - a)))/2;
 }
 
 
@@ -68,14 +68,15 @@ Foam::blendingMethods::hyperbolic::hyperbolic
     ),
     transitionAlphaScale_
     (
-        readParameter("transitionAlphaScale", dict, {0, NaN}, false)
+        readParameter("transitionAlphaScale", dict, {0, vGreat}, false)
     )
 {
     if
     (
         canBeContinuous(0)
      && canBeContinuous(1)
-     && minContinuousAlpha_[0] + minContinuousAlpha_[1] < 1 - rootSmall
+     && minContinuousAlpha_[0].value + minContinuousAlpha_[1].value
+      < 1 - rootSmall
     )
     {
         FatalErrorInFunction
@@ -98,7 +99,7 @@ Foam::blendingMethods::hyperbolic::~hyperbolic()
 
 bool Foam::blendingMethods::hyperbolic::canBeContinuous(const label index) const
 {
-    return isParameter(minContinuousAlpha_[index]);
+    return minContinuousAlpha_[index].valid;
 }
 
 
@@ -107,7 +108,8 @@ bool Foam::blendingMethods::hyperbolic::canSegregate() const
     return
         canBeContinuous(0)
      && canBeContinuous(1)
-     && minContinuousAlpha_[0] + minContinuousAlpha_[1] > 1 + rootSmall;
+     && minContinuousAlpha_[0].value + minContinuousAlpha_[1].value
+      > 1 + rootSmall;
 }
 
 

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,10 +65,10 @@ class hyperbolic
     // Private Data
 
         //- Minimum fraction of phases which can be considered continuous
-        Pair<scalar> minContinuousAlpha_;
+        Pair<blendingParameter> minContinuousAlpha_;
 
         //- Width of the transition
-        const scalar transitionAlphaScale_;
+        const blendingParameter transitionAlphaScale_;
 
 
 protected:

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,8 +95,8 @@ Foam::blendingMethods::linear::linear
 
         if
         (
-            isParameter(minFullyContinuousAlpha_[i])
+            minFullyContinuousAlpha_[i].valid
-         != isParameter(minPartlyContinuousAlpha_[i])
+         != minPartlyContinuousAlpha_[i].valid
         )
         {
             FatalErrorInFunction
@@ -109,7 +109,8 @@ Foam::blendingMethods::linear::linear
         (
             (
                 canBeContinuous(i)
-             && minFullyContinuousAlpha_[i] <= minPartlyContinuousAlpha_[i]
+             && minFullyContinuousAlpha_[i].value
+             <= minPartlyContinuousAlpha_[i].value
             )
         )
         {
@@ -126,11 +127,13 @@ Foam::blendingMethods::linear::linear
      && canBeContinuous(1)
      && (
             (
-                minFullyContinuousAlpha_[0] + minPartlyContinuousAlpha_[1]
+                minFullyContinuousAlpha_[0].value
+              + minPartlyContinuousAlpha_[1].value
               < 1 - rootSmall
             )
          || (
-                minFullyContinuousAlpha_[1] + minPartlyContinuousAlpha_[0]
+                minFullyContinuousAlpha_[1].value
+              + minPartlyContinuousAlpha_[0].value
               < 1 - rootSmall
             )
         )
@@ -157,7 +160,7 @@ Foam::blendingMethods::linear::~linear()
 
 bool Foam::blendingMethods::linear::canBeContinuous(const label index) const
 {
-    return isParameter(minFullyContinuousAlpha_[index]);
+    return minFullyContinuousAlpha_[index].valid;
 }
 
 
@@ -167,11 +170,16 @@ bool Foam::blendingMethods::linear::canSegregate() const
         canBeContinuous(0)
      && canBeContinuous(1)
      && (
-            minFullyContinuousAlpha_[0] + minPartlyContinuousAlpha_[1]
+            (
-          > 1 + rootSmall
+                minFullyContinuousAlpha_[0].value
-         ||
+              + minPartlyContinuousAlpha_[1].value
-            minFullyContinuousAlpha_[1] + minPartlyContinuousAlpha_[0]
+              > 1 + rootSmall
-          > 1 + rootSmall
+            )
+         || (
+                minFullyContinuousAlpha_[1].value
+              + minPartlyContinuousAlpha_[0].value
+              > 1 + rootSmall
+            )
         );
 }
 

applications/solvers/modules/multiphaseEuler/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,10 +64,10 @@ class linear
     // Private Data
 
         //- Minimum fraction of phases which can be considered fully continuous
-        Pair<scalar> minFullyContinuousAlpha_;
+        Pair<blendingParameter> minFullyContinuousAlpha_;
 
         //- Minimum fraction of phases which can be considered partly continuous
-        Pair<scalar> minPartlyContinuousAlpha_;
+        Pair<blendingParameter> minPartlyContinuousAlpha_;
 
 
 protected: