diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
index f714be3..9439fc4 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
@@ -64,27 +64,42 @@ Foam::BreakupModel<CloudType>::BreakupModel
 (
     const dictionary& dict,
     CloudType& owner,
-    const word& type
+    const word& type,
+    const Switch solveOscillationEq = false
 )
 :
     SubModelBase<CloudType>(owner, dict, typeName, type),
-    solveOscillationEq_(this->coeffDict().lookup("solveOscillationEq")),
+    solveOscillationEq_(solveOscillationEq),
     y0_(0.0),
     yDot0_(0.0),
-    TABComega_(0.0),
-    TABCmu_(0.0),
-    TABWeCrit_(0.0)
+    TABComega_(8.0),
+    TABCmu_(10.0),
+    TABWeCrit_(12.0)
 {
+    // Give the option to solve the oscillation equation, even if not required 
+    // by the breakup model
+    if (!solveOscillationEq_)
+    {
+        this->coeffDict().readIfPresent
+        (
+            "solveOscillationEq",
+            solveOscillationEq_
+        );
+    }
+
+    // Check for non-default TAB coefficients
     if (solveOscillationEq_)
     {
         const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
-        y0_ = TABcoeffsDict.template lookupOrDefault<scalar>("y0", 0.0);
-        yDot0_ = TABcoeffsDict.template lookupOrDefault<scalar>("yDot0", 0.0);
-        TABComega_ =
-            TABcoeffsDict.template lookupOrDefault<scalar>("Comega", 8.0);
-        TABCmu_ = TABcoeffsDict.template lookupOrDefault<scalar>("Cmu", 10.0);
-        TABWeCrit_ =
-            TABcoeffsDict.template lookupOrDefault<scalar>("WeCrit", 12.0);
+        Switch defaultCoeffs = TABcoeffsDict.lookup("defaultCoeffs");
+        if (!defaultCoeffs)
+        {
+            TABcoeffsDict.lookup("y0") >> y0_;
+            TABcoeffsDict.lookup("yDot0") >> yDot0_;
+            TABcoeffsDict.lookup("Comega") >> TABComega_;
+            TABcoeffsDict.lookup("Cmu") >> TABCmu_;
+            TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
+        }
     }
 }
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
index e0672ca..09cbe36 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
@@ -97,7 +97,8 @@ public:
         (
             const dictionary& dict,
             CloudType& owner,
-            const word& type
+            const word& type,
+            const Switch solveOscillationEq
         );
 
         //- Construct copy
diff --git a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
index fb48040..7f08e47 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
@@ -34,7 +34,7 @@ Foam::ETAB<CloudType>::ETAB
     CloudType& owner
 )
 :
-    BreakupModel<CloudType>(dict, owner, typeName),
+    BreakupModel<CloudType>(dict, owner, typeName, true),
     Cmu_(10.0),
     Comega_(8.0),
     k1_(0.2),
@@ -55,15 +55,6 @@ Foam::ETAB<CloudType>::ETAB
 
     scalar k21 = k2_/k1_;
     AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
-
-    if (!BreakupModel<CloudType>::solveOscillationEq_)
-    {
-        Info<< "Warning: solveOscillationEq is set to "
-            << BreakupModel<CloudType>::solveOscillationEq_ << nl
-            << " Setting it to true in order for the ETAB model to work."
-            << endl;
-        BreakupModel<CloudType>::solveOscillationEq_ = true;
-    }
 }
 
 
diff --git a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
index 0fe325f..53f416f 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/TAB/TAB.C
@@ -34,7 +34,7 @@ Foam::TAB<CloudType>::TAB
     CloudType& owner
 )
 :
-    BreakupModel<CloudType>(dict, owner, typeName),
+    BreakupModel<CloudType>(dict, owner, typeName, true),
     Cmu_(BreakupModel<CloudType>::TABCmu_),
     Comega_(BreakupModel<CloudType>::TABComega_),
     WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
@@ -52,16 +52,6 @@ Foam::TAB<CloudType>::TAB
             (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
     }
 
-    if (!BreakupModel<CloudType>::solveOscillationEq_)
-    {
-        WarningIn("Foam::TAB<CloudType>::TAB(const dictionary&, CloudType&)")
-            << "solveOscillationEq is set to "
-            << BreakupModel<CloudType>::solveOscillationEq_ << nl
-            << " Setting it to true in order for the TAB model to work."
-            << endl;
-        BreakupModel<CloudType>::solveOscillationEq_ = true;
-    }
-
     if (SMDCalcMethod_ == "method1")
     {
         SMDMethod_ = method1;

From c0e01d022e1923d080872c812408e5f78165ebe2 Mon Sep 17 00:00:00 2001
From: Armin Wehrfritz <armin.wehrfritz@aalto.fi>
Date: Wed, 4 Sep 2013 21:32:39 +0300
Subject: [PATCH] Fix inconsistent ETAB coefficient reading

---
 src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C |    9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
index 7f08e47..24a33e1 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C
@@ -35,21 +35,18 @@ Foam::ETAB<CloudType>::ETAB
 )
 :
     BreakupModel<CloudType>(dict, owner, typeName, true),
-    Cmu_(10.0),
-    Comega_(8.0),
+    Cmu_(BreakupModel<CloudType>::TABCmu_),
+    Comega_(BreakupModel<CloudType>::TABComega_),
     k1_(0.2),
     k2_(0.2),
-    WeCrit_(12.0),
+    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
     WeTransition_(100.0),
     AWe_(0.0)
 {
     if (!this->defaultCoeffs(true))
     {
-        this->coeffDict().lookup("Cmu") >> Cmu_;
-        this->coeffDict().lookup("Comega") >> Comega_;
         this->coeffDict().lookup("k1") >> k1_;
         this->coeffDict().lookup("k2") >> k2_;
-        this->coeffDict().lookup("WeCrit") >> WeCrit_;
         this->coeffDict().lookup("WeTransition") >> WeTransition_;
     }
 
-- 
1.7.10.4

From 48872ed9df56504a8b3a60eff752ff758cc75049 Mon Sep 17 00:00:00 2001
From: Armin Wehrfritz <armin.wehrfritz@aalto.fi>
Date: Wed, 4 Sep 2013 22:26:42 +0300
Subject: [PATCH] Enhance TAB coefficient reading and add
 enableOscillationEq()

---
 .../BreakupModel/BreakupModel/BreakupModel.C       |   54 ++++++++++++++++----
 .../BreakupModel/BreakupModel/BreakupModel.H       |   10 ++++
 2 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
index 9439fc4..f5fa020 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C
@@ -25,6 +25,30 @@ License
 
 #include "BreakupModel.H"
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::BreakupModel<CloudType>::readTABCoeffs(const bool printMsg)
+{
+    const dictionary TABcoeffsDict(this->dict().subDict("TABCoeffs"));
+    bool def = TABcoeffsDict.lookupOrDefault<bool>("defaultCoeffs", false);
+    if (!def)
+    {
+        TABcoeffsDict.lookup("y0") >> y0_;
+        TABcoeffsDict.lookup("yDot0") >> yDot0_;
+        TABcoeffsDict.lookup("Comega") >> TABComega_;
+        TABcoeffsDict.lookup("Cmu") >> TABCmu_;
+        TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
+    }
+    if (printMsg && def)
+    {
+        Info<< incrIndent;
+        Info<< indent << "Employing default TAB coefficients" << endl;
+        Info<< decrIndent;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -90,16 +114,7 @@ Foam::BreakupModel<CloudType>::BreakupModel
     // Check for non-default TAB coefficients
     if (solveOscillationEq_)
     {
-        const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
-        Switch defaultCoeffs = TABcoeffsDict.lookup("defaultCoeffs");
-        if (!defaultCoeffs)
-        {
-            TABcoeffsDict.lookup("y0") >> y0_;
-            TABcoeffsDict.lookup("yDot0") >> yDot0_;
-            TABcoeffsDict.lookup("Comega") >> TABComega_;
-            TABcoeffsDict.lookup("Cmu") >> TABCmu_;
-            TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
-        }
+        readTABCoeffs(true);
     }
 }
 
@@ -114,6 +129,25 @@ Foam::BreakupModel<CloudType>::~BreakupModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
+void Foam::BreakupModel<CloudType>::enableOscillationEq(const bool printMsg)
+{
+    if (!solveOscillationEq_)
+    {
+        solveOscillationEq_ = true;
+
+        if (printMsg)
+        {
+            Info<< incrIndent;
+            Info<< indent << "Enable oscillation equation solution" << endl;
+            Info<< decrIndent;
+        }
+
+        readTABCoeffs(printMsg);
+    }
+}
+
+
+template<class CloudType>
 bool Foam::BreakupModel<CloudType>::update
 (
     const scalar dt,
diff --git a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
index 09cbe36..ac54b73 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H
@@ -55,6 +55,12 @@ class BreakupModel
     public SubModelBase<CloudType>
 {
 
+    // Private Member Functions
+
+        //- Read TAB coefficients
+        void readTABCoeffs(const bool printMsg);
+
+
 protected:
 
     // Protected data
@@ -161,6 +167,10 @@ public:
 
     // Member Functions
 
+        //- Enable the solution of the oscillation equation and check that the 
+        //  TAB coefficients are set correctly
+        void enableOscillationEq(const bool printMsg);
+
         //- Update the parcel properties and return true if a child parcel
         //  should be added
         virtual bool update
-- 
1.7.10.4

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::BreakupModel

Description
    Templated break-up model class

SourceFiles
    BreakupModel.C
    BreakupModelNew.C

\*---------------------------------------------------------------------------*/

#ifndef BreakupModel_H
#define BreakupModel_H

#include "IOdictionary.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                       Class BreakupModel Declaration
\*---------------------------------------------------------------------------*/

template<class CloudType>
class BreakupModel
:
    public SubModelBase<CloudType>
{

    // Private Member Functions

        //- Read TAB coefficients
        void readTABCoeffs(const bool printMsg);


protected:

    // Protected data

        Switch solveOscillationEq_;

        scalar y0_;
        scalar yDot0_;
        scalar TABComega_;
        scalar TABCmu_;
        scalar TABWeCrit_;


public:

    //- Runtime type information
    TypeName("breakupModel");

    //- Declare runtime constructor selection table
    declareRunTimeSelectionTable
    (
        autoPtr,
        BreakupModel,
        dictionary,
        (
            const dictionary& dict,
            CloudType& owner
        ),
        (dict, owner)
    );


    // Constructors

        //- Construct null from owner
        BreakupModel(CloudType& owner);

        //- Construct from dictionary
        BreakupModel
        (
            const dictionary& dict,
            CloudType& owner,
            const word& type,
            const Switch solveOscillationEq
        );

        //- Construct copy
        BreakupModel(const BreakupModel<CloudType>& bum);

        //- Construct and return a clone
        virtual autoPtr<BreakupModel<CloudType> > clone() const
        {
            return autoPtr<BreakupModel<CloudType> >
            (
                new BreakupModel<CloudType>(*this)
            );
        }


    //- Destructor
    virtual ~BreakupModel();


    //- Selector
    static autoPtr<BreakupModel<CloudType> > New
    (
        const dictionary& dict,
        CloudType& owner
    );


    // Access

        inline const Switch& solveOscillationEq() const
        {
            return solveOscillationEq_;
        }

        inline const scalar& y0() const
        {
            return y0_;
        }

        inline const scalar& yDot0() const
        {
            return yDot0_;
        }

        inline const scalar& TABComega() const
        {
            return TABComega_;
        }

        inline const scalar& TABCmu() const
        {
            return TABCmu_;
        }

        inline const scalar& TABWeCrit() const
        {
            return TABWeCrit_;
        }


    // Member Functions

        //- Enable the solution of the oscillation equation and check that the 
        //  TAB coefficients are set correctly
        void enableOscillationEq(const bool printMsg);

        //- Update the parcel properties and return true if a child parcel
        //  should be added
        virtual bool update
        (
            const scalar dt,
            const vector& g,
            scalar& d,
            scalar& tc,
            scalar& ms,
            scalar& nParticle,
            scalar& KHindex,
            scalar& y,
            scalar& yDot,
            const scalar d0,
            const scalar rho,
            const scalar mu,
            const scalar sigma,
            const vector& U,
            const scalar rhoc,
            const scalar muc,
            const vector& Urel,
            const scalar Urmag,
            const scalar tMom,
            scalar& dChild,
            scalar& massChild
        );
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#define makeBreakupModel(CloudType)                                           \
                                                                              \
    typedef CloudType::sprayCloudType sprayCloudType;                         \
    defineNamedTemplateTypeNameAndDebug                                       \
    (                                                                         \
        BreakupModel<sprayCloudType>,                                         \
        0                                                                     \
    );                                                                        \
    defineTemplateRunTimeSelectionTable                                       \
    (                                                                         \
        BreakupModel<sprayCloudType>,                                         \
        dictionary                                                            \
    );


#define makeBreakupModelType(SS, CloudType)                                   \
                                                                              \
    typedef CloudType::sprayCloudType sprayCloudType;                         \
    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
                                                                              \
    BreakupModel<sprayCloudType>::                                            \
        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
            add##SS##CloudType##sprayCloudType##ConstructorToTable_;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "BreakupModel.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "BreakupModel.H"

// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //

template<class CloudType>
void Foam::BreakupModel<CloudType>::readTABCoeffs(const bool printMsg)
{
    const dictionary TABcoeffsDict(this->dict().subDict("TABCoeffs"));
    bool def = TABcoeffsDict.lookupOrDefault<bool>("defaultCoeffs", false);
    if (!def)
    {
        TABcoeffsDict.lookup("y0") >> y0_;
        TABcoeffsDict.lookup("yDot0") >> yDot0_;
        TABcoeffsDict.lookup("Comega") >> TABComega_;
        TABcoeffsDict.lookup("Cmu") >> TABCmu_;
        TABcoeffsDict.lookup("WeCrit") >> TABWeCrit_;
    }
    if (printMsg && def)
    {
        Info<< incrIndent;
        Info<< indent << "Employing default TAB coefficients" << endl;
        Info<< decrIndent;
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    CloudType& owner
)
:
    SubModelBase<CloudType>(owner),
    solveOscillationEq_(false),
    y0_(0.0),
    yDot0_(0.0),
    TABComega_(0.0),
    TABCmu_(0.0),
    TABWeCrit_(0.0)
{}


template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    const BreakupModel<CloudType>& bum
)
:
    SubModelBase<CloudType>(bum),
    solveOscillationEq_(bum.solveOscillationEq_),
    y0_(bum.y0_),
    yDot0_(bum.yDot0_),
    TABComega_(bum.TABComega_),
    TABCmu_(bum.TABCmu_),
    TABWeCrit_(bum.TABWeCrit_)
{}


template<class CloudType>
Foam::BreakupModel<CloudType>::BreakupModel
(
    const dictionary& dict,
    CloudType& owner,
    const word& type,
    const Switch solveOscillationEq = false
)
:
    SubModelBase<CloudType>(owner, dict, typeName, type),
    solveOscillationEq_(solveOscillationEq),
    y0_(0.0),
    yDot0_(0.0),
    TABComega_(8.0),
    TABCmu_(10.0),
    TABWeCrit_(12.0)
{
    // Give the option to solve the oscillation equation, even if not required 
    // by the breakup model
    if (!solveOscillationEq_)
    {
        this->coeffDict().readIfPresent
        (
            "solveOscillationEq",
            solveOscillationEq_
        );
    }

    // Check for non-default TAB coefficients
    if (solveOscillationEq_)
    {
        readTABCoeffs(true);
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::BreakupModel<CloudType>::~BreakupModel()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
void Foam::BreakupModel<CloudType>::enableOscillationEq(const bool printMsg)
{
    if (!solveOscillationEq_)
    {
        solveOscillationEq_ = true;

        if (printMsg)
        {
            Info<< incrIndent;
            Info<< indent << "Enable oscillation equation solution" << endl;
            Info<< decrIndent;
        }

        readTABCoeffs(printMsg);
    }
}


template<class CloudType>
bool Foam::BreakupModel<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    notImplemented
    (
        "bool Foam::BreakupModel<CloudType>::update"
        "("
            "const scalar, "
            "const vector&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "scalar&, "
            "const scalar, "
            "const scalar, "
            "const scalar, "
            "const scalar, "
            "const vector&, "
            "const scalar, "
            "const scalar, "
            "const vector&, "
            "const scalar, "
            "const scalar, "
            "scalar&, "
            "scalar&"
        ");"
    );

    return false;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "BreakupModelNew.C"

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "ETAB.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::ETAB<CloudType>::ETAB
(
    const dictionary& dict,
    CloudType& owner
)
:
    BreakupModel<CloudType>(dict, owner, typeName, true),
    Cmu_(BreakupModel<CloudType>::TABCmu_),
    Comega_(BreakupModel<CloudType>::TABComega_),
    k1_(0.2),
    k2_(0.2),
    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
    WeTransition_(100.0),
    AWe_(0.0)
{
    if (!this->defaultCoeffs(true))
    {
        this->coeffDict().lookup("k1") >> k1_;
        this->coeffDict().lookup("k2") >> k2_;
        this->coeffDict().lookup("WeTransition") >> WeTransition_;
    }

    scalar k21 = k2_/k1_;
    AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
}


template<class CloudType>
Foam::ETAB<CloudType>::ETAB(const ETAB<CloudType>& bum)
:
    BreakupModel<CloudType>(bum),
    Cmu_(bum.Cmu_),
    Comega_(bum.Comega_),
    k1_(bum.k1_),
    k2_(bum.k2_),
    WeCrit_(bum.WeCrit_),
    WeTransition_(bum.WeTransition_),
    AWe_(bum.AWe_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::ETAB<CloudType>::~ETAB()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
bool Foam::ETAB<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    scalar r  = 0.5*d;
    scalar r2 = r*r;
    scalar r3 = r*r2;

    scalar semiMass = nParticle*pow3(d);

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar romega = 1.0/omega;

        scalar We = rhoc*sqr(Urmag)*r/sigma;
        scalar Wetmp = We/WeCrit_;

        scalar y1 = y - Wetmp;
        scalar y2 = yDot*romega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a + Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = constant::mathematical::twoPi - phit;
            }

            scalar tb = 0;

            if (mag(y) < 1.0)
            {
                scalar theta = acos((1.0 - Wetmp)/a);

                if (theta < phi)
                {
                    if (constant::mathematical::twoPi - theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += constant::mathematical::twoPi;
                }
                tb = (theta - phi)*romega;

                // breakup occurs
                if (dt > tb)
                {
                    y = 1.0;
                    yDot = -a*omega*sin(omega*tb + phi);
                }
            }

            // update droplet size
            if (dt > tb)
            {
                scalar sqrtWe = AWe_*pow4(We) + 1.0;
                scalar Kbr = k1_*omega*sqrtWe;

                if (We > WeTransition_)
                {
                    sqrtWe = sqrt(We);
                    Kbr =k2_*omega*sqrtWe;
                }

                scalar rWetmp = 1.0/Wetmp;
                scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
                scalar dtbu = romega*acos(cosdtbu);
                scalar decay = exp(-Kbr*dtbu);

                scalar rNew = decay*r;
                if (rNew < r)
                {
                    d = 2.0*rNew;
                    y = 0.0;
                    yDot = 0.0;
                }
            }
        }
    }
    else
    {
        // reset droplet distortion parameters
        y = 0;
        yDot = 0;
    }

    // update the nParticle count to conserve mass
    nParticle = semiMass/pow3(d);

    // Do not add child parcel
    return false;
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "TAB.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::TAB<CloudType>::TAB
(
    const dictionary& dict,
    CloudType& owner
)
:
    BreakupModel<CloudType>(dict, owner, typeName, true),
    Cmu_(BreakupModel<CloudType>::TABCmu_),
    Comega_(BreakupModel<CloudType>::TABComega_),
    WeCrit_(BreakupModel<CloudType>::TABWeCrit_),
    SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
{
    // calculate the inverse function of the Rossin-Rammler Distribution
    const scalar xx0 = 12.0;
    const scalar rrd100 =
        1.0/(1.0 - exp(-xx0)*(1.0 + xx0 + sqr(xx0)/2.0 + pow3(xx0)/6.0));

    forAll(rrd_, n)
    {
        scalar xx = 0.12*(n + 1);
        rrd_[n] =
            (1.0 - exp(-xx)*(1.0 + xx + sqr(xx)/2.0 + pow3(xx)/6.0))*rrd100;
    }

    if (SMDCalcMethod_ == "method1")
    {
        SMDMethod_ = method1;
    }
    else if (SMDCalcMethod_ == "method2")
    {
        SMDMethod_ = method2;
    }
    else
    {
        SMDMethod_ = method2;
        WarningIn("Foam::TAB<CloudType>::TAB(const dictionary&, CloudType&)")
            << "Unknown SMDCalculationMethod. Valid options are "
            << "(method1 | method2). Using method2" << endl;
    }
}


template<class CloudType>
Foam::TAB<CloudType>::TAB(const TAB<CloudType>& bum)
:
    BreakupModel<CloudType>(bum),
    Cmu_(bum.Cmu_),
    Comega_(bum.Comega_),
    WeCrit_(bum.WeCrit_),
    SMDCalcMethod_(bum.SMDCalcMethod_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CloudType>
Foam::TAB<CloudType>::~TAB()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class CloudType>
bool Foam::TAB<CloudType>::update
(
    const scalar dt,
    const vector& g,
    scalar& d,
    scalar& tc,
    scalar& ms,
    scalar& nParticle,
    scalar& KHindex,
    scalar& y,
    scalar& yDot,
    const scalar d0,
    const scalar rho,
    const scalar mu,
    const scalar sigma,
    const vector& U,
    const scalar rhoc,
    const scalar muc,
    const vector& Urel,
    const scalar Urmag,
    const scalar tMom,
    scalar& dChild,
    scalar& massChild
)
{
    cachedRandom& rndGen = this->owner().rndGen();

    scalar r = 0.5*d;
    scalar r2 = r*r;
    scalar r3 = r*r2;

    scalar semiMass = nParticle*pow3(d);

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar We = rhoc*sqr(Urmag)*r/sigma;
        scalar Wetmp = We/WeCrit_;

        scalar y1 = y - Wetmp;
        scalar y2 = yDot/omega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a+Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = constant::mathematical::twoPi - phit;
            }

            scalar tb = 0;

            if (mag(y) < 1.0)
            {
                scalar coste = 1.0;
                if ((Wetmp - a < -1) && (yDot < 0))
                {
                    coste = -1.0;
                }

                scalar theta = acos((coste-Wetmp)/a);

                if (theta < phi)
                {
                    if (constant::mathematical::twoPi - theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += constant::mathematical::twoPi;
                }
                tb = (theta-phi)/omega;

                // breakup occurs
                if (dt > tb)
                {
                    y = 1.0;
                    yDot = -a*omega*sin(omega*tb + phi);
                }

            }

            // update droplet size
            if (dt > tb)
            {
                scalar rs =
                    r/(1.0 + (4.0/3.0)*sqr(y) + rho*r3/(8*sigma)*sqr(yDot));

                label n = 0;
                scalar rNew = 0.0;
                switch (SMDMethod_)
                {
                    case method1:
                    {
                        #include "TABSMDCalcMethod1.H"
                        break;
                    }
                    case method2:
                    {
                        #include "TABSMDCalcMethod2.H"
                        break;
                    }
                }

                if (rNew < r)
                {
                    d = 2*rNew;
                    y = 0;
                    yDot = 0;
                }
            }
        }
    }
    else
    {
        // reset droplet distortion parameters
        y = 0;
        yDot = 0;
    }

    // update the nParticle count to conserve mass
    nParticle = semiMass/pow3(d);

    // Do not add child parcel
    return false;
}


// ************************************************************************* //