LagrangianCMCFoam AMC integration algorithm update

solvers_post/BetaIntegrator/BetaIntegrator.H

@@ -69,15 +69,15 @@ class BetaIntegrator
 {
     // Private data
 
+        //- Quadrature points
+        const scalarField &etaSpace_;
+
         //- Beta pdf numerators
         scalarField pdfNum_;
 
         //- Beta pdf denominator
         scalar pdfDen_;
 
-        //- Quadrature points
-        const scalarField &etaSpace_;
-
         //- Description of data_
         BetaFunction bf_;
 

solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C

@@ -51,7 +51,13 @@ Contact
 
 #include "LinearMesh.H"
 
-#include "../FlameStructure/FlameStructure.H"
+#include "FlameStructure.H"
+
+#include "AMC.H"
+
+#include "BetaFunction.H"
+#include "BetaGrid.H"
+#include "BetaIntegrator.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

solvers_post/LagrangianCMCFoam/Make/files

@@ -1,5 +1,20 @@
 ../FlameStructure/FlameStructure.C
 ../FlameStructure/FlameStructureIO.C
+
+../AMC/AMC.C
+
+../LinearInterpolator1D/LinearInterpolator1D.C
+../LinearInterpolator1D/LinearInterpolator1DIO.C
+
+../BetaFunction/BetaFunction.C
+../BetaFunction/BetaFunctionIO.C
+
+../BetaGrid/BetaGrid.C
+../BetaGrid/BetaGridIO.C
+
+../BetaIntegrator/BetaIntegrator.C
+../BetaIntegrator/BetaIntegratorIO.C
+
 LagrangianCMCFoam.C
 
 EXE = $(FOAM_USER_APPBIN)/LagrangianCMCFoam

solvers_post/LagrangianCMCFoam/Make/options

@@ -1,6 +1,12 @@
 DEV_PATH=../../libs
 
 EXE_INC = \
+    -I../FlameStructure \
+    -I../LinearInterpolator1D \
+    -I../AMC \
+    -I../BetaFunction \
+    -I../BetaGrid \
+    -I../BetaIntegrator \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \

solvers_post/LagrangianCMCFoam/readCMCProperties.H

@@ -95,6 +95,9 @@ scalarField pdf(etamax+1,0.0);
 scalarField f(etamax+1, 0.0);
 Info<<"Set eta space grid"<<tab<<uni_eta<<endl;
 
+
+CMC::BetaGrid bg(CMCdict.subDict("betaGrid"));
+
 //Uniform eta space grid
 if(uni_eta == true)
 {

solvers_post/LagrangianCMCFoam/setSDRCMC.H

@@ -1,36 +1,34 @@
-//get conditional SDR
-        scalar C2(0); 
-        scalar Airside(0);
-        scalar Fuelside(0);
+// Calculate conditional SDR
 
-        forAll(rho , celli) 
+        scalarField amc(CMC::AMC(bg.etaSpace()));
+
+        forAll(rho, celli)
         {
-                Airside = Peta[0][celli]*deta[0]; 
-                Fuelside = Peta[etamax][celli]*deta[etamax]; 
-                if( (1-Airside) < 0.000001 || (1-Fuelside) < 0.000001 )//if there exist delta function at eta = 0 or 1
-                {
-                        for(label j = 0 ; j<= etamax ; j++)
-                        {
-                                CSDRCELL[j][celli] = 0;
-                        }
-                }
-                else
-                {
-                        scalar C1(0);
-			for(label j=0 ; j<=etamax ;j++)
-			{
-				f[j] = AMC(etaValue[j]);
-				pdf[j] = Peta[j][celli];
-			}
-			C1 = integration(deltaftn[celli], MFcut, Neta, pdf, f);
-                        C2 = SDR[celli] / C1; //Eqn (13) from SH Kim C&F paper 120:75-90(20)
+            scalar mfi = mf[celli];
+            scalar mfVari = mfVar[celli];
 
-                        for(label j = 0 ; j <= etamax ; j++)
-                        {
-                                scalar eta = etaValue[j];
-                                CSDRCELL[j][celli] = C2 * AMC(eta);
-                        }
+            CMC::BetaFunction bf (mfi, mfVari);
+            CMC::BetaIntegrator bi (bf, bg);
+
+            //if there exist delta function at eta = 0 or 1
+            if( bf.delta0() || bf.delta1() )
+            {
+                for(label j = 0 ; j<= etamax ; j++)
+                {
+                    CSDRCELL[j][celli] = 0;
                 }
+            }
+            else
+            {
+                //Eqn (13) from SH Kim C&F paper 120:75-90(20)
+                scalar C1 = bi.betaIntegrate(amc);
+                scalar C2 = SDR[celli] / C1;
+
+                for(label j = 0 ; j <= etamax ; j++)
+                {
+                    CSDRCELL[j][celli] = C2 * CMC::AMC(etaValue[j]);
+                }
+            }
         }
 
 //get density weighted averaged value of conditional SDR for each eta and flame-group 

solvers_post/LinearInterpolator1D/LinearInterpolator1D.C

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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 "LinearInterpolator1D.H"
+
+#include "error.H"
+#include "Switch.H"
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+bool Foam::LinearInterpolator1D::_checkSorted(const scalarField& grid)
+{
+    Switch sorted = true;
+
+    for (label i = 1; i < srcGrid_.size(); i++)
+    {
+        if (grid[i-1] >= grid[i])
+        {
+            sorted = false;
+            break;
+        }
+    }
+
+    return sorted;
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LinearInterpolator1D::LinearInterpolator1D()
+{}
+
+
+Foam::LinearInterpolator1D::LinearInterpolator1D(const scalarField& srcGrid, const scalarField& dstGrid)
+:
+    srcGrid_(srcGrid),
+    dstGrid_(dstGrid),
+    idxHigh_(dstGrid.size()),
+    coef_(dstGrid.size())
+{
+    // Check dstGrid and srcGrid are intervals
+    if (!(dstGrid_.first() < dstGrid_.last() && srcGrid_.first() < srcGrid_.last()))
+    {
+        FatalErrorInFunction
+            << "Invalid grids"
+            << abort(FatalError);
+    }
+
+    // Check dstGrid lies within srcGrid
+    if (dstGrid_.first() < srcGrid_.first() || srcGrid_.last() < dstGrid_.last())
+    {
+        FatalErrorInFunction
+            << "Extrapolation attempted"
+            << abort(FatalError);
+    }
+
+    // Check srcGrid and dstGrid are sorted
+    if (!(_checkSorted(dstGrid_) && _checkSorted(srcGrid_)))
+    {
+        FatalErrorInFunction
+            << "Unsorted grids"
+            << abort(FatalError);
+    }
+
+    forAll (dstGrid_, j)
+    {
+        for (label i = 1; i < srcGrid_.size(); i++)
+        {
+            if (srcGrid_[i] >= dstGrid_[j])
+            {
+                idxHigh_[j] = i;
+                coef_[j] = (srcGrid_[i] - dstGrid_[j]) / (srcGrid_[i] - srcGrid_[i-1]);
+                break;
+            }
+        }
+    }
+}
+
+
+Foam::LinearInterpolator1D::LinearInterpolator1D(const LinearInterpolator1D&)
+{}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::LinearInterpolator1D::~LinearInterpolator1D()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+void Foam::LinearInterpolator1D::interpolate(const scalarField& src, scalarField& dst)
+{
+    forAll (dst, i)
+    {
+        dst[i] = coef_[i]*src[idxHigh_[i]-1] + (1.0-coef_[i])*src[idxHigh_[i]];
+    }
+
+    return;
+}
+
+
+// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
+
+void Foam::LinearInterpolator1D::operator=(const LinearInterpolator1D& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+}
+
+// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

solvers_post/LinearInterpolator1D/LinearInterpolator1D.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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::LinearInterpolator1D
+
+Description
+
+SourceFiles
+    LinearInterpolator1DI.H
+    LinearInterpolator1D.C
+    LinearInterpolator1DIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LinearInterpolator1D_H
+#define LinearInterpolator1D_H
+
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class LinearInterpolator1D;
+Istream& operator>>(Istream&, LinearInterpolator1D&);
+Ostream& operator<<(Ostream&, const LinearInterpolator1D&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class LinearInterpolator1D Declaration
+\*---------------------------------------------------------------------------*/
+
+class LinearInterpolator1D
+{
+    // Private data
+
+        //- Source grids
+        scalarField srcGrid_;
+
+        //- Destination grids
+        scalarField dstGrid_;
+
+        //- Source grid to refer for each dst grid
+        labelList idxHigh_;
+
+        //- Coefficients for dst grid points
+        scalarField coef_;
+
+
+    // Private Member Functions
+
+        //- Check grid sorted
+        bool _checkSorted(const scalarField& grid);
+
+        //- Disallow default bitwise copy construct
+        LinearInterpolator1D(const LinearInterpolator1D&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const LinearInterpolator1D&);
+
+
+public:
+
+
+    // Constructors
+
+        //- Construct null
+        LinearInterpolator1D();
+
+        //- Construct from components
+        LinearInterpolator1D(const scalarField& srcGrid, const scalarField& dstGrid);
+
+        //- Construct from Istream
+        LinearInterpolator1D(Istream&);
+
+
+
+    //- Destructor
+    ~LinearInterpolator1D();
+
+
+    // Member Functions
+
+        // Access
+
+        // Check
+
+        // Edit
+
+        // Write
+
+        //- Evaluate interpolation
+        void interpolate(const scalarField& src, scalarField& dst);
+
+
+    // Member Operators
+
+
+
+    // Friend Functions
+
+    // Friend Operators
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, LinearInterpolator1D&);
+        friend Ostream& operator<<(Ostream&, const LinearInterpolator1D&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "LinearInterpolator1DI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

solvers_post/LinearInterpolator1D/LinearInterpolator1DI.H

@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

solvers_post/LinearInterpolator1D/LinearInterpolator1DIO.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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 "LinearInterpolator1D.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LinearInterpolator1D::LinearInterpolator1D(Istream& is)
+{
+    // Check state of Istream
+    is.check("Foam::LinearInterpolator1D::LinearInterpolator1D(Foam::Istream&)");
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, LinearInterpolator1D&)
+{
+    // Check state of Istream
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::LinearInterpolator1D&)"
+    );
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const LinearInterpolator1D&)
+{
+    // Check state of Ostream
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
+        "const Foam::LinearInterpolator1D&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //