From a8044e1f535ea910bfe2dc32b5278513d6e515ce Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Tue, 23 Jan 2018 21:35:04 +0900
Subject: [PATCH] remove CMC methods and pure virtual methods added to
 chemistryModels

---
 .../basicChemistryModel/basicChemistryModel.H |  34 +-
 .../chemistryModel/chemistryModel.C           | 353 +-----------------
 .../chemistryModel/chemistryModel.H           |  29 --
 3 files changed, 8 insertions(+), 408 deletions(-)

diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index 2d636fd89..553fc1a01 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -178,45 +178,23 @@ public:
 
                 //- Return the mass fraction of the *specieName
                 //- 28.Nov.2017 Jinwoo Park
-                virtual tmp<volScalarField> Yspecie(const word *specieName) const = 0;
+                virtual tmp<volScalarField> Yspecie(const word *specieName) {}
 
                 //- Return the mass fraction of the *specieName
                 //- 28.Nov.2017 Jinwoo Park
-                virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) const = 0;
+                virtual tmp<volScalarField> MoleFracSpecie(const word *specieName) {}
 
                 //- Return the epsilon
-                virtual tmp<volScalarField> epsilon() const = 0;
+                virtual tmp<volScalarField> epsilon() {}
 
                 //- Return the k
-                virtual tmp<volScalarField> k() const = 0;
+                virtual tmp<volScalarField> k() {}
 
                 //- Return the U field
-                virtual tmp<volVectorField> U() const = 0;
+                virtual tmp<volVectorField> U() {}
 
                 //- Return the turbulent kinematic diffusivity of chemical species (Le=1 is assumed)
-                virtual tmp<volScalarField> nut() const = 0;
-
-                virtual scalar solveCMCchem
-                (
-                    scalar deltaT, scalar& rho, scalar& T,
-                    scalar& p, scalarField& Qi, scalar& Qh,
-                    scalarField& RRCMC, scalar& ChemDeltaT
-                ) = 0; //CMC
-
-                virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p) = 0; //CMC
-                virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p) = 0; //CMC
-                virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p) = 0; //CMC
-                virtual scalar calculateShCMC(scalarField& RRCMC, label i) = 0; //CMC
-            virtual void correction
-            (
-                scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
-                scalarField& Gab, scalarField& Gat, scalarField& Gbt,
-                scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
-            ) = 0;
-
-                //- Calculate the reaction rate (EDM or finite-rate/EDM)
-                //- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
-                virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B) = 0;
+                virtual tmp<volScalarField> nut() {}
 
 };
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 6d99a3c01..6807df96e 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -27,7 +27,6 @@ License
 #include "reactingMixture.H"
 #include "UniformField.H"
 #include "extrapolatedCalculatedFvPatchFields.H"
-#include "wallFvPatch.H"                                             //jinwoo, implemented according to the of231 version of simpleCoalCombustionfoam from karam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -50,12 +49,11 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
         dynamic_cast<const reactingMixture<ThermoType>&>
             (this->thermo()).speciesData()
     ),
+
     nSpecie_(Y_.size()),
     nReaction_(reactions_.size()),
     Treact_(CompType::template lookupOrDefault<scalar>("Treact", 0.0)),
-    RR_(nSpecie_),
-    edmRR_(reactions_.size()),
-    kineticRR_(reactions_.size())
+    RR_(nSpecie_)
 {
     // create the fields for the chemistry sources
     forAll(RR_, fieldi)
@@ -79,45 +77,6 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
         );
     }
 
-    forAll(edmRR_, reactioni)
-    {
-        edmRR_.set
-        (
-            reactioni,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "edmRR." + reactions_[reactioni].name(),
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh,
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-            )
-        );
-
-        kineticRR_.set
-        (
-            reactioni,
-            new volScalarField
-            (
-                IOobject
-                (
-                    "kineticRR." + reactions_[reactioni].name(),
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh,
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-            )
-        );
-    }
-
     Info<< "chemistryModel: Number of species = " << nSpecie_
         << " and reactions = " << nReaction_ << endl;
 }
@@ -1161,313 +1120,5 @@ void Foam::chemistryModel<CompType, ThermoType>::solve
     NotImplemented;
 }
 
-template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solveCMCchem
-(
-    scalar deltaT, scalar& rho, scalar& T, scalar& p, scalarField& Qi, scalar& Qh, scalarField& RRCMC, scalar& ChemDeltaT
-)
-{
-    scalar deltaTMin = GREAT;
-    //if(!this->chemistry_)
-    //{
-    //    return deltaTMin;
-    //}
-
-    scalar rhoi = rho;
-    scalar Ti = T;
-    scalar pi = p;
-
-    scalarField c(nSpecie_);
-    scalarField c0(nSpecie_);
-
-    for(label i=0 ; i<nSpecie_ ; i++)
-    {
-        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
-        c0[i] = c[i];
-    }
-
-    scalar timeLeft = deltaT;
-
-    while(timeLeft > SMALL)
-    {
-        scalar dt = timeLeft;
-        this->solve(c, Ti, pi, dt, ChemDeltaT);
-        timeLeft -= dt;
-    }
-    deltaTMin = min(ChemDeltaT, deltaTMin);
-
-    for(label i=0 ; i<nSpecie_; i++)
-    {
-        RRCMC[i] = (c[i]-c0[i])*specieThermo_[i].W()/deltaT;
-    }
-    return deltaTMin;
-}
-
-template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateTCMC
-(
-    scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p
-)
-{
-    scalar rhoi = rho;
-    scalar Ti = Told;
-    scalar pi = p;
-
-    scalarField c(nSpecie_);
-
-    for (label i=0; i<nSpecie_; i++)
-    {
-        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
-    }
-
-    scalar cTot = 0.0;
-
-    // update the temperature
-    cTot = sum(c);
-    ThermoType mixture(0.0*specieThermo_[0]);
-    for (label i=0; i<nSpecie_; i++)
-    {
-        mixture += (c[i]/cTot)*specieThermo_[i];
-    }
-    scalar TCMC = mixture.THa(Qh, pi,  Ti);
-    //scalar TCMCnew = mixture.THs(Qh, Ti);
-
-    return TCMC;//
-}
-
-template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateRHOCMC
-(
-    scalarField& Qi, scalar& T, scalar& p
-)
-{
-    scalar Ti = T; //conditional temperature [K]
-    scalar pi = p; //pressure [Pa]
-
-    scalar MeanMW(0); //mean molecular weight [kg/kmol]
-    scalar InvMeanMW(0); //inverse of MeanMW [kmol/kg]
-    scalar RHOCMC(0); //conditional density [kg/kmol]
-
-    for (label i=0; i<nSpecie_; i++) //get mean molecular weight respect to Qi
-    {
-        InvMeanMW += (Qi[i]/specieThermo_[i].W());
-    }
-    MeanMW = 1.0/InvMeanMW ;
-
-    RHOCMC = pi/(8314.4621/MeanMW*Ti); //universial gas constant 8314.51[J/kmol K] = [kg.m2/kmol.K.sec2]
-
-    return RHOCMC;//
-}
-
-template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateHCMC
-(
-    scalarField& Qi, scalar& T, scalar& rho, scalar& p
-)
-{
-    scalar rhoi = rho; //conditional density [kg/kmol]
-    scalar pi = p;
-    scalar Ti = T; //conditional temperature [K]
-    scalar HCMC(0); //conditional total enthalpy [J/kg]
-    scalarField c(nSpecie_);
-
-    for (label i=0; i<nSpecie_; i++)
-    {
-        c[i] = rhoi*Qi[i]/specieThermo_[i].W();
-    }
-
-    scalar cTot = 0.0;
-
-    // update the temperature
-    cTot = sum(c);
-    ThermoType mixture(0.0*specieThermo_[0]);
-    for (label i=0; i<nSpecie_; i++)
-    {
-        mixture += (c[i]/cTot)*specieThermo_[i];
-    }
-    HCMC = mixture.Ha(pi, Ti);
-
-    return HCMC;
-}
-
-template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::calculateShCMC
-(
-    scalarField& RRCMC, label i
-)
-{
-    scalar tSh(0);
-    if (this->chemistry_)
-    {
-        scalar hi = specieThermo_[i].Hc();
-        //Info<<"hi = "<<hi<<endl;
-        //Info<<"RRCMC[i] = "<<RRCMC[i]<<endl;
-        tSh = hi*RRCMC[i];
-        //Info<<"tSh = "<<tSh<<endl;
-    }
-    return tSh;
-}
-
-template<class CompType, class ThermoType>
-void Foam::chemistryModel<CompType, ThermoType>::correction
-(
-        scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
-        scalarField& Gab, scalarField& Gat, scalarField& Gbt,
-        scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
-)
-{
-        for(label i=0 ; i<3 ; i++)
-        {
-                b_[i] = b[i];
-                Ta_[i] = Ta[i];
-                Wa_[i] = Wa[i];
-                Wb_[i] = Wb[i];
-
-                Gab_[i] = Gab[i];
-                Gat_[i] = Gat[i];
-                Gbt_[i] = Gbt[i];
-                Gaa_[i] = Gaa[i];
-                Gbb_[i] = Gbb[i];
-        }
-        Gtt_ = Gtt;
-        rhorho_ = rhorho;
-}
-
-/*
-template<class CompType, class ThermoType>
-void Foam::chemistryModel<CompType, ThermoType>::specieMoleFrac
-(
-    char specieName
-)
-{
-}
-*/
-
-template<class CompType, class ThermoType>
-void Foam::chemistryModel<CompType, ThermoType>::calculateEDM
-(
-    bool finiteRate, const scalar A, const scalar B
-)
-{
-    if (!this->chemistry_)
-    {
-        Info<<"Reaction rates are not calculated!!"<<endl;
-        Info<<"Please check constant/chemistryProperties"<<endl;
-        return;
-    }
-
-    const volScalarField& rho = this->thermo().rho();
-    const volScalarField& k = this->db().objectRegistry::lookupObject<volScalarField>("k");
-    const volScalarField& epsilon = this->db().objectRegistry::lookupObject<volScalarField>("epsilon");
-
-    const scalarField& T = this->thermo().T();
-    const scalarField& p = this->thermo().p();
-
-    scalar EDMrate(0.0);
-
-    forAll(rho, celli)
-    {
-        scalarField RRedm(nSpecie_, 0.0);
-        scalarField RRkinetic(nSpecie_, 0.0);
-
-        forAll(reactions_, reactioni)
-        {
-            const Reaction<ThermoType>& R = reactions_[reactioni];
-
-            scalar Ymin(1.0);
-            forAll(R.lhs(), lhs_speciei)
-            {
-                const label speciei = R.lhs()[lhs_speciei].index;
-                const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
-                if(stoiCoeff != 0)
-                {
-                    Ymin = min(Ymin, Y_[speciei][celli]/stoiCoeff/specieThermo_[speciei].W());
-                }
-            }
-
-            scalar YPmin(0.0);
-            scalar Num(0.0);
-            scalar Den(0.0);
-            forAll(R.rhs(), rhs_speciei)
-            {
-                const label speciei = R.rhs()[rhs_speciei].index;
-                const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
-                if(stoiCoeff != 0)
-                {
-                    Num += Y_[speciei][celli];
-                    Den += stoiCoeff*specieThermo_[speciei].W();
-                }
-            }
-            YPmin = Num / Den;
-
-            Ymin = min(Ymin, B*YPmin);
-
-            EDMrate = A*rho[celli]*epsilon[celli]/k[celli]*Ymin;
-
-            scalar omegai = 0.0;
-            if(finiteRate == true)
-            {
-                scalar pf, cf, pr, cr;
-                label lRef, rRef;
-
-                scalarField c(nSpecie_, 0.0);
-                for(label i=0 ; i<nSpecie_ ; i++)
-                {
-                    const scalar Yi = Y_[i][celli];
-                    c[i] = rho[celli]*Yi/specieThermo_[i].W();
-                }
-                omegai = omega
-                (
-                    R, c, T[celli], p[celli], pf, cf, lRef, pr, cr, rRef
-                );
-            }
-
-            forAll(R.lhs(), lhs_speciei)
-            {
-                const label speciei = R.lhs()[lhs_speciei].index;
-                const scalar stoiCoeff = R.lhs()[lhs_speciei].stoichCoeff;
-
-                RRedm[speciei] -= stoiCoeff*specieThermo_[speciei].W()*EDMrate;
-                RRkinetic[speciei] -= stoiCoeff*specieThermo_[speciei].W()*omegai;
-
-                if(lhs_speciei == 0)
-                {
-                    //for post-processing
-                    edmRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*EDMrate;
-                    kineticRR_[reactioni][celli] = stoiCoeff*specieThermo_[speciei].W()*omegai;
-                }
-            }
-            forAll(R.rhs(), rhs_speciei)
-            {
-                const label speciei = R.rhs()[rhs_speciei].index;
-                const scalar stoiCoeff = R.rhs()[rhs_speciei].stoichCoeff;
-
-                RRedm[speciei] += stoiCoeff*specieThermo_[speciei].W()*EDMrate;
-                RRkinetic[speciei] += stoiCoeff*specieThermo_[speciei].W()*omegai;
-
-            }
-        }
-
-        for (label i=0; i<nSpecie_; i++)
-        {
-            if(finiteRate == true)
-            {
-                if(RRedm[i] == 0.0 || RRkinetic[i] == 0.0)
-                {
-                    RR_[i][celli] = 0.0;
-                }
-                else
-                {
-                    RR_[i][celli] = (RRedm[i]*RRkinetic[i])/(RRedm[i]+RRkinetic[i]);
-                }
-            }
-            else
-            {
-                RR_[i][celli] = RRedm[i];
-            }
-        }
-    }
-}
-
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 1228e75e3..d5343adea 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -103,13 +103,6 @@ protected:
         //- List of reaction rate per specie [kg/m3/s]
         PtrList<DimensionedField<scalar, volMesh>> RR_;
 
-        scalarField b_, Ta_, Wa_, Wb_, Gab_, Gat_, Gbt_, Gaa_, Gbb_;
-
-        scalar Gtt_, rhorho_;
-
-        PtrList<DimensionedField<scalar, volMesh>> edmRR_;
-
-        PtrList<DimensionedField<scalar, volMesh>> kineticRR_;
 
     // Protected Member Functions
 
@@ -229,24 +222,6 @@ public:
             //  and return the characteristic time
             virtual scalar solve(const scalarField& deltaT);
 
-            virtual scalar solveCMCchem
-            (
-                scalar deltaT, scalar& rho, scalar& T,
-                scalar& p, scalarField& Qi, scalar& Qh,
-                scalarField& RRCMC, scalar& ChemDeltaT
-            );
-
-            virtual scalar calculateTCMC(scalar& Qh, scalarField& Qi, scalar& Told, scalar& rho, scalar& p);
-            virtual scalar calculateRHOCMC(scalarField& Qi, scalar& T, scalar& p);
-            virtual scalar calculateHCMC(scalarField& Qi, scalar& T, scalar& rho, scalar& p);
-            virtual scalar calculateShCMC(scalarField& RRCMC, label i);
-            virtual void correction
-            (
-                scalarField& b, scalarField& Ta, scalarField& Wa, scalarField& Wb,
-                scalarField& Gab, scalarField& Gat, scalarField& Gbt,
-                scalarField& Gaa, scalarField& Gbb, scalar& Gtt, scalar& rhorho
-            );
-
             //- Return the chemical time scale
             virtual tmp<volScalarField> tc() const;
 
@@ -280,10 +255,6 @@ public:
             //- 28.Nov.2017 Jinwoo Park
             virtual tmp<volScalarField> nut() const;
 
-            //- Calculate the reaction rates (EDM or finite-rate EDM)
-            //- 11.Apr.2017 Karam Han (updated by Jinwoo Park)
-            virtual void calculateEDM(bool finiteRate, const scalar A, const scalar B);
-
         // ODE functions (overriding abstract functions in ODE.H)
 
             //- Number of ODE's to solve