recurrence solver for modal decomposition of flow fleids

applications/solvers/rctfSpeciesTransport/CEq.H

@@ -0,0 +1,19 @@
+
+volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + alphat);
+
+CEqn =
+(
+    fvm::ddt(C)
+    + fvm::div(phiRec, C)
+    - fvm::laplacian(alphaEff, C)
+    ==
+    fvOptions(C)
+);
+
+CEqn.relax(relaxCoeff);
+
+fvOptions.constrain(CEqn);
+
+CEqn.solve();
+
+fvOptions.correct(C);

applications/solvers/rctfSpeciesTransport/Make/files

@@ -0,0 +1,3 @@
+rctfSpeciesTransport.C
+
+EXE=$(CFDEM_APP_DIR)/rctfSpeciesTransport

applications/solvers/rctfSpeciesTransport/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/rctfSpeciesTransport/continuityErrCalc.H

@@ -0,0 +1,17 @@
+// calculate the continuity error according to phiRec
+
+{
+    volScalarField contErr(fvc::div(phiRec));
+
+    scalar sumLocalContErr = runTime.deltaTValue()*
+        mag(contErr)().weightedAverage(mesh.V()).value();
+
+    scalar globalContErr = runTime.deltaTValue()*
+        contErr.weightedAverage(mesh.V()).value();
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}

applications/solvers/rctfSpeciesTransport/createFields.H

@@ -0,0 +1,281 @@
+//creating the fields according to the recurrence dictionary
+
+IOdictionary recProperties_
+(
+    IOobject
+    (
+        "recProperties0",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+List<wordList> fieldsDict_(recProperties_.lookup("fieldsPairs"));
+
+wordList fieldNames(fieldsDict_.size());
+
+for(int i = 0; i < fieldsDict_.size(); i++)
+{
+
+    fieldNames[i]= fieldsDict_[i][0];
+}
+
+Info<< "\n list of the fields: \n" << fieldNames << endl;
+
+//reading coherent velocity field name
+label k = findIndex(fieldNames,"coh_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent velocity\n" << abort(FatalError);
+}
+const word Ucoh_pair = fieldsDict_[k][1];
+
+    volVectorField UcohRec
+    (
+        IOobject
+        (
+            Ucoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading incoherent velocity field name
+k = findIndex(fieldNames,"inc_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the incoherent velocity\n" << abort(FatalError);
+}
+const word Uinc_pair = fieldsDict_[k][1];
+
+    volVectorField UincRec
+    (
+        IOobject
+        (
+            Uinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading coherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_coh");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGScoh_pair = fieldsDict_[k][1];
+
+    volScalarField kcohRec
+    (
+        IOobject
+        (
+            kSGScoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+//reading incoherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_inc");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGSinc_pair = fieldsDict_[k][1];
+
+    volScalarField kincRec
+    (
+        IOobject
+        (
+            kSGSinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+// calculated fields
+Info<< "\nCreating cell volume field\n" << endl;
+
+volScalarField delta
+(
+    IOobject
+    (
+        "delta",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+   mesh,
+   dimensionedScalar("delta", dimLength, 0.0)
+);
+
+delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0);
+delta.write();
+
+
+volVectorField URec
+(
+    IOobject
+    (
+        "URec",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+Info<< "\nCreating turb kinetic energy field\n" << endl;
+
+volScalarField kRec
+(
+    IOobject
+    (
+        "kRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    kcohRec+kincRec
+);
+
+// check if there is any negative values
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+
+kRec.write();
+
+Info<< "\nCreating turb viscosity field\n" << endl;
+volScalarField nutRec
+(
+    IOobject
+    (
+        "nutRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    sqrt(kRec)*delta*0.094
+);
+
+nutRec.write();
+
+
+Info<< "Calculating face flux field phiRec\n" << endl;
+surfaceScalarField phiRec
+ (
+     IOobject
+     (
+        "phiRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(URec) & mesh.Sf()
+ );
+
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+  
+    dimensionedScalar Sc("Sc", dimless, laminarTransport);
+    dimensionedScalar Sct("Sct", dimless, laminarTransport);
+  
+volScalarField alphat
+(
+    IOobject
+    (
+        "alphat",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    nutRec/Sct
+);
+
+ // create the scalar field
+ Info<< "Creating scalar transport field\n" << endl;
+
+    volScalarField C
+    (
+        IOobject
+        (
+            "C",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    fvScalarMatrix CEqn(C, dimless*dimVolume/(dimTime));
+
+    scalar relaxCoeff(0.0);
+
+ Info<< "reading clockProperties\n" << endl;
+
+    IOdictionary clockProperties
+    (
+            IOobject
+            (
+                "clockProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+    );
+
+    autoPtr<clockModel> myClock
+    (
+            clockModel::New
+            (
+                clockProperties,
+                mesh.time()
+            )
+    );

applications/solvers/rctfSpeciesTransport/createRecBase.H

@@ -0,0 +1,66 @@
+// check which recProperties dicts are present, read them in and construct a PtrList of recBases
+// names for dicts can be "recProperties" or "recPropertiesN" where N in {0, 1, ...} 
+
+#include "error.H"
+
+word dictName = "recProperties";
+wordList recPropertiesList(0);
+PtrList <recBase> recBases(0);
+label maxDictNumber = 100;
+
+{
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictName,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictName);
+    }
+}
+
+
+for (label counter = 0; counter < maxDictNumber; counter++)
+{
+    word dictNameIter = dictName + Foam::name(counter);
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictNameIter,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictNameIter);
+    }
+}
+
+if (recPropertiesList.size() == 0)
+{
+    FatalError << "no recProperties dicts found" << endl;
+}
+
+else
+{
+    Info << "found " << recPropertiesList.size() << " dicts with names " << recPropertiesList << endl;
+}
+
+
+for (label counter = 0; counter < recPropertiesList.size(); counter++)
+{
+    recBases.append( new recBase(mesh, recPropertiesList[counter]));
+}

applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger, Sanaz Abbasi
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Turbulent Transport Recurrence Solver for modal decomposition
+
+Description
+    Solves a transport equation for a passive scalar on a single-phase solution
+    for a solver based on recurrence statistics
+
+Rules
+	Solution data to compute the recurrence statistics from, needs to 
+        reside in $CASE_ROOT/dataBase(0...N)
+    Time step data in the first dataBase needs to be evenly spaced in time
+    A list of indices for the corresponding incoherent fields to coherent ones
+        should be provided.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "clockModel.H"
+
+#include "objectRegistry.H"
+#include "VectorSpace.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    scalar cumulativeContErr = 0;
+  
+    //create recBases according to a list of recProperties
+    #include "createRecBase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+    label recTimeIndex(0);
+    label currTimeIndex(0);
+
+    scalar recTimeStep_=recBases[0].recM().recTimeStep();
+    labelPairList incPairTimeIndex_(0);
+
+    IFstream pairFile("incIndexPairList");
+    pairFile >> incPairTimeIndex_;
+
+    while (runTime.run())
+    {
+
+        myClock().start(1,"Global");
+        runTime++;
+
+        myClock().start(11,"Total");
+        
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        myClock().start(2,"fieldUpdate");
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	    Info << "Updating fields at run time " << runTime.timeOutputValue()
+	        << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recBases[0].updateRecFields();
+            #include "readFields.H"
+
+            recTimeIndex++;
+        }
+
+        myClock().stop("fieldUpdate");
+
+        #include "continuityErrCalc.H"
+
+        myClock().start(3,"speciesEqn");
+        #include "CEq.H"
+        myClock().stop("speciesEqn");
+
+        myClock().stop("Total");
+
+        runTime.write();
+        
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+        myClock().stop("Global");
+
+    }
+
+    
+    myClock().evalPar();
+    myClock().normHist();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rctfSpeciesTransport/readFields.H

@@ -0,0 +1,40 @@
+
+currTimeIndex = recBases[0].recM().currentTimeIndex();
+
+Info << "current Time Index = " << currTimeIndex << endl;
+
+recBases[0].recM().exportVolVectorField(Ucoh_pair,UcohRec);
+recBases[0].recM().exportVolScalarField(kSGScoh_pair,kcohRec);
+
+
+label incTimeIndex = incPairTimeIndex_[currTimeIndex][1];
+
+Info << "  incoherent pair Time Index = " << incTimeIndex  << endl;
+
+UincRec = recBases[1].recM().exportVolVectorField(Uinc_pair,incTimeIndex);
+kincRec = recBases[1].recM().exportVolScalarField(kSGSinc_pair,incTimeIndex);
+
+kRec = kcohRec+kincRec;
+
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+URec = UcohRec + UincRec;
+phiRec = linearInterpolate(URec) & mesh.Sf();
+
+nutRec = sqrt(kRec)*delta*0.094;
+
+alphat = nutRec/Sct;
+alphat.correctBoundaryConditions();
+

etc/solver-list.txt

@@ -12,3 +12,4 @@ cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
 cfdemSolverRhoPimpleChem/dir
 cfdemSolverMultiphase/dir
+rctfSpeciesTransport/dir