Bulk.hpp

Go to the documentation of this file.

 
#ifndef __BULK_HEADER__
#define __BULK_HEADER__
 
// Standard header files
#include <iostream>
#include <vector>
 
// OpenCAEPoro header files
#include "DenseMat.hpp"
#include "FlowUnit.hpp"
#include "Grid.hpp"
#include "LinearSystem.hpp"
#include "Mixture.hpp"
#include "MixtureBO.hpp"
#include "MixtureComp.hpp"
#include "MixtureThermal.hpp"
#include "OCPConst.hpp"
#include "OCPStructure.hpp"
#include "ParamReservoir.hpp"
#include "Rock.hpp"
 
using namespace std;
 
// General Error Yype
 
const int BULK_SUCCESS                   = 0;
const int BULK_NEGATIVE_PRESSURE         = -1;
const int BULK_NEGATIVE_TEMPERATURE      = -2;
const int BULK_NEGATIVE_COMPONENTS_MOLES = -3;
const int BULK_OUTRANGED_VOLUME_ERROR    = -4;
const int BULK_OUTRANGED_CFL             = -5;
 
//  Note: This class includes reference depth and pressure at it, depth of contacts
//  between phases, and capillary pressure at phase contact surfaces.
class ParamEQUIL
{
    friend class Bulk;
 
private:
    OCP_DBL  Dref; 
    OCP_DBL  Pref; 
    OCP_DBL  DOWC; 
    OCP_DBL  DGOC; 
    OCP_DBL  PcOW; 
    OCP_DBL  PcGO; 
    OCPTable PBVD; 
};
 
class HeatLoss
{
public:
    void     InputParam(const HLoss& loss);
    void     Setup(const OCP_USI& nb);
    OCP_BOOL IfHeatLoss() const { return ifHLoss; }
    void     CalHeatLoss(const vector<USI>&     location,
                         const vector<OCP_DBL>& T,
                         const vector<OCP_DBL>& lT,
                         const vector<OCP_DBL>& initT,
                         const OCP_DBL&         t,
                         const OCP_DBL&         dt);
    void     ResetToLastTimeStep();
    void     UpdateLastTimeStep();
 
public:
    OCP_BOOL ifHLoss{OCP_FALSE}; 
    OCP_DBL  obC;                
    OCP_DBL  obK;                
    OCP_DBL  ubC;                
    OCP_DBL  ubK;                
 
    OCP_USI         numBulk; 
    OCP_DBL         obD;     
    OCP_DBL         ubD;     
    vector<OCP_DBL> I;       
    vector<OCP_DBL> p;       
    vector<OCP_DBL> pT;      
 
    vector<OCP_DBL> lI;  
    vector<OCP_DBL> lp;  
    vector<OCP_DBL> lpT; 
};
 
//  Note: Bulk contains main physical infomation of active grids. It describes the
//  actual geometric domain for simulating. Variables are stored bulk by bulk, and then
//  phase by phase, then component by component. The bulks are ordered in the alphabetic
//  order, i.e. the X-axis indices first, followed by the Y-axis and Z-axis indices.
//  Operations on each bulk are also defined here.
class Bulk
{
    friend class BulkConn;
    friend class Well;
    friend class Out4RPT;
    friend class Out4VTK;
 
    // temp
    friend class Reservoir;
    friend class IsothermalMethod;
    friend class IsoT_IMPEC;
    friend class IsoT_FIM;
    friend class IsoT_FIMn;
    friend class IsoT_AIMc;
    friend class T_FIM;
 
    // Input Param and Setup
 
public:
    void InputParam(const ParamReservoir& rs_param);
    void InputParamBLKOIL(const ParamReservoir& rs_param);
    void InputParamCOMPS(const ParamReservoir& rs_param);
    void InputParamTHERMAL(const ParamReservoir& rs_param);
    void InputSatFunc(const ParamReservoir& rs_param);
    void InputRockFunc(const ParamReservoir& rs_param);
    void InputRockFuncT(const ParamReservoir& rs_param);
 
    void SetupIsoT(const Grid& myGrid);
    void SetupT(const Grid& myGrid);
    void SetupOptionalFeatures(const Grid& myGrid, OptionalFeatures& optFeatures);
 
    // General Variables
 
public:
    OCP_USI GetBulkNum() const { return numBulk; }
    USI GetPhaseNum() const { return numPhase; }
    USI GetComNum() const { return numCom; }
 
protected:
    OCP_USI numBulk;  
    USI     numPhase; 
    USI     numCom;   
    USI     numComH;  
 
    // Initial Properties
public:
    void InitPTSw(const USI& tabrow);
 
protected:
    vector<OCPTable>
        initZi_Tab; 
    vector<OCPTable> initT_Tab; 
    vector<OCP_DBL>  initT;     
    ParamEQUIL       EQUIL;     
    OCP_DBL          rsTemp;    
    vector<OCP_DBL>  thconp;    
 
    // Region
 
public:
    void AllocateRegion(const Grid& myGrid);
    void SetupBulkType(const Grid& myGrid);
    const vector<Mixture*>& GetMixture() const { return flashCal; }
    void OutMixtureIters() const { flashCal[0]->OutMixtureIters(); }
 
protected:
    USI              NTPVT;    
    USI              PVTmodeB; 
    vector<USI>      PVTNUM;   
    vector<Mixture*> flashCal; 
 
    USI               NTSFUN;  
    USI               SATmode; 
    vector<USI>       SATNUM;  
    vector<FlowUnit*> flow;    
    vector<vector<OCP_DBL>>
        satcm; 
 
    USI           NTROCC;  
    vector<USI>   ROCKNUM; 
    vector<Rock*> rock;    
 
    vector<USI> bType;     
    vector<USI> bLocation; 
    HeatLoss    hLoss;     
 
    // Basic PVT Model Information
 
public:
    OCP_BOOL IfUseEoS() const { return ifUseEoS; };
 
protected:
    OCP_BOOL ifBlackOil{OCP_FALSE}; 
    OCP_BOOL ifComps{OCP_FALSE};    
    OCP_BOOL ifThermal{OCP_FALSE};  
    OCP_BOOL ifUseEoS{OCP_FALSE};   
    OCP_BOOL oil{OCP_FALSE};        
    OCP_BOOL gas{OCP_FALSE};        
    OCP_BOOL water{OCP_FALSE};      
    OCP_BOOL disGas{OCP_FALSE}; 
 
    // Basic Grid and Basic Rock Information
 
public:
    void AllocateGridRockIsoT(const Grid& myGrid);
    void AllocateGridRockT(const Grid& myGrid);
 
protected:
    vector<OCP_DBL> dx;    
    vector<OCP_DBL> dy;    
    vector<OCP_DBL> dz;    
    vector<OCP_DBL> v;     
    vector<OCP_DBL> depth; 
 
    vector<OCP_DBL> ntg;      
    vector<OCP_DBL> poroInit; 
    vector<OCP_DBL> poro;     
    vector<OCP_DBL> rockVp;   
    vector<OCP_DBL> rockKx;   
    vector<OCP_DBL> rockKy;   
    vector<OCP_DBL> rockKz;   
    vector<OCP_DBL> thconr;   
    vector<OCP_DBL> vr;       
    vector<OCP_DBL> Hr;       
 
    // Last time step
    vector<OCP_DBL> lporo;   
    vector<OCP_DBL> lrockVp; 
    vector<OCP_DBL> lvr;     
    vector<OCP_DBL> lHr;     
 
    // Derivatives
    vector<OCP_DBL> poroP; 
    vector<OCP_DBL> poroT; 
    vector<OCP_DBL> vrP;   
    vector<OCP_DBL> vrT;   
    vector<OCP_DBL> HrT;   
 
    // Last time step
    vector<OCP_DBL> lporoP; 
    vector<OCP_DBL> lporoT; 
    vector<OCP_DBL> lvrP;   
    vector<OCP_DBL> lvrT;   
    vector<OCP_DBL> lHrT;   
 
    // Basic Fluid Information
 
public:
    OCP_DBL CalFPR() const;
    OCP_DBL CalFTR() const;
    OCP_DBL GetP(const OCP_USI& n) const { return P[n]; }
    OCP_DBL GetSOIL(const OCP_USI& n) const
    {
        return S[n * numPhase + phase2Index[OIL]];
    }
    OCP_DBL GetSGAS(const OCP_USI& n) const
    {
        return S[n * numPhase + phase2Index[GAS]];
    }
    OCP_DBL GetSWAT(const OCP_USI& n) const
    {
        return S[n * numPhase + phase2Index[WATER]];
    }
 
protected:
    vector<USI> phase2Index;     
                                 // Note: For example, `Oil' is at the i-th location.
    vector<USI>      phaseNum;   
    vector<OCP_DBL>  Nt;         
    vector<OCP_DBL>  Ni;         
    vector<OCP_DBL>  vf;         
    vector<OCP_DBL>  T;          
    vector<OCP_DBL>  P;          
    vector<OCP_DBL>  Pb;         
    vector<OCP_DBL>  Pj;         
    vector<OCP_DBL>  Pc;         
    vector<OCP_BOOL> phaseExist; 
    vector<OCP_DBL>  S;          
    vector<OCP_DBL>  vj;         
    vector<OCP_DBL>  nj;         
    vector<OCP_DBL>  xij;        
    vector<OCP_DBL>  rho;        
    vector<OCP_DBL>  xi;         
    vector<OCP_DBL>  mu;         
    vector<OCP_DBL>  kr;         
    vector<OCP_DBL>  Uf;         
    vector<OCP_DBL>  H;          
    vector<OCP_DBL>  kt;         
 
    // Last time step
    vector<USI>      lphaseNum;   
    vector<OCP_DBL>  lNt;         
    vector<OCP_DBL>  lNi;         
    vector<OCP_DBL>  lvf;         
    vector<OCP_DBL>  lT;          
    vector<OCP_DBL>  lP;          
    vector<OCP_DBL>  lPj;         
    vector<OCP_DBL>  lPc;         
    vector<OCP_BOOL> lphaseExist; 
    vector<OCP_DBL>  lS;          
    vector<OCP_DBL>  lvj;         
    vector<OCP_DBL>  lnj;         
    vector<OCP_DBL>  lxij;        
    vector<OCP_DBL>  lrho;        
    vector<OCP_DBL>  lxi;         
    vector<OCP_DBL>  lmu;         
    vector<OCP_DBL>  lkr;         
    vector<OCP_DBL>  lUf;         
    vector<OCP_DBL>  lH;          
    vector<OCP_DBL>  lkt;         
 
    // Derivatives
    vector<OCP_DBL> vfP;     
    vector<OCP_DBL> vfT;     
    vector<OCP_DBL> vfi;     
    vector<OCP_DBL> rhoP;    
    vector<OCP_DBL> rhoT;    
    vector<OCP_DBL> rhox;    
    vector<OCP_DBL> xiP;     
    vector<OCP_DBL> xiT;     
    vector<OCP_DBL> xix;     
    vector<OCP_DBL> muP;     
    vector<OCP_DBL> muT;     
    vector<OCP_DBL> mux;     
    vector<OCP_DBL> dPcj_dS; 
    vector<OCP_DBL> dKr_dS;  
    vector<OCP_DBL> UfP;     
    vector<OCP_DBL> UfT;     
    vector<OCP_DBL> Ufi;     
    vector<OCP_DBL> HT;      
    vector<OCP_DBL> Hx;      
    vector<OCP_DBL> ktP;     
    vector<OCP_DBL> ktT;     
    vector<OCP_DBL> ktS;     
 
    // Last time step
    vector<OCP_DBL> lvfP;     
    vector<OCP_DBL> lvfT;     
    vector<OCP_DBL> lvfi;     
    vector<OCP_DBL> lrhoP;    
    vector<OCP_DBL> lrhoT;    
    vector<OCP_DBL> lrhox;    
    vector<OCP_DBL> lxiP;     
    vector<OCP_DBL> lxiT;     
    vector<OCP_DBL> lxix;     
    vector<OCP_DBL> lmuP;     
    vector<OCP_DBL> lmuT;     
    vector<OCP_DBL> lmux;     
    vector<OCP_DBL> ldPcj_dS; 
    vector<OCP_DBL> ldKr_dS;  
    vector<OCP_DBL> lUfP;     
    vector<OCP_DBL> lUfT;     
    vector<OCP_DBL> lUfi;     
    vector<OCP_DBL> lHT;      
    vector<OCP_DBL> lHx;      
    vector<OCP_DBL> lktP;     
    vector<OCP_DBL> lktT;     
    vector<OCP_DBL> lktS;     
 
    // Newton Iteration Information
 
public:
    OCP_DBL CalNRdSmax(OCP_USI& index);
    OCP_DBL GetNRdPmax() const { return NRdPmax; };
    OCP_DBL GetNRdNmax() const { return NRdNmax; };
 
protected:
    vector<OCP_DBL> dSNR;  
    vector<OCP_DBL> dSNRP; 
    vector<OCP_DBL> dNNR;  
    vector<OCP_DBL> dPNR;  
    vector<OCP_DBL> dTNR;  
 
    OCP_DBL maxNRdSSP;       
    OCP_USI index_maxNRdSSP; 
    OCP_DBL NRdPmax;         
    OCP_DBL NRdTmax;         
    OCP_DBL NRdNmax;         
    OCP_DBL NRdSmax;         
 
    vector<OCP_DBL> NRstep;     
    vector<USI>     NRphaseNum; 
 
    // Important Indicator Variable and Check
 
public:
    OCP_INT CheckP() const;
    OCP_INT CheckT() const;
    OCP_INT CheckNi();
    OCP_INT CheckVe(const OCP_DBL& Vlim) const;
    OCP_INT CheckCFL(const OCP_DBL& cflLim) const;
 
    void CalMaxChange();
 
    OCP_DBL GetdPmax() const { return dPmax; }
    OCP_DBL GetdTmax() const { return dTmax; }
    OCP_DBL GetdNmax() const { return dNmax; }
    OCP_DBL GetdSmax() const { return dSmax; }
    OCP_DBL GeteVmax() const { return eVmax; }
    OCP_DBL GetMaxCFL() const { return maxCFL; }
 
protected:
    OCP_DBL dPmax; 
    OCP_DBL dTmax; 
    OCP_DBL dSmax; 
    OCP_DBL dNmax; 
    OCP_DBL eVmax; 
 
    mutable vector<OCP_DBL> cfl;    
    mutable OCP_DBL         maxCFL; 
 
    // Error
 
public:
    void AllocateError();
 
protected:
    vector<OCP_DBL> ePEC; 
    mutable OCPRes  res;  
 
protected:
    // Method-Specified Variable
 
    USI              maxLendSdP;   
    vector<USI>      bRowSizedSdP; 
    vector<OCP_DBL>  dSec_dPri;    
    vector<OCP_BOOL> pSderExist;   
    vector<USI>      pVnumCom;     
    vector<OCP_DBL>  res_n;        
    vector<OCP_DBL>  resPc;        
 
    // Last time step
    vector<USI>      lbRowSizedSdP; 
    vector<OCP_DBL>  ldSec_dPri;    
    vector<OCP_BOOL> lpSderExist;   
    vector<USI>      lpVnumCom;     
    vector<OCP_DBL>  lres_n;        
    vector<OCP_DBL>  lresPc;        
 
public:
    void ShowFIMBulk(const OCP_BOOL& flag = OCP_FALSE) const;
    void AddWellBulkId(const OCP_USI& n) { wellBulkId.push_back(n); }
 
protected:
    vector<OCP_USI> wellBulkId; 
    class BulkTypeAIM
    {
    public:
        void Setup(const OCP_USI& nb) { indicator.resize(nb, -1); }
        void Init() { fill(indicator.begin(), indicator.end(), -1); }
        void SetBulkType(const OCP_USI& n, const OCP_INT& flag) { indicator[n] = flag; }
        OCP_BOOL IfFIMbulk(const OCP_USI& n) const { return indicator[n] > 0; }
        OCP_BOOL IfIMPECbulk(const OCP_USI& n) const { return indicator[n] < 0; }
 
    protected:
        vector<OCP_INT> indicator; 
    } bulkTypeAIM;
    vector<OCP_DBL> xijNR; 
};
 
#endif /* end if __BULK_HEADER__ */
 
/*----------------------------------------------------------------------------*/
/*  Brief Change History of This File                                         */
/*----------------------------------------------------------------------------*/
/*  Author              Date             Actions                              */
/*----------------------------------------------------------------------------*/
/*  Shizhe Li           Oct/01/2021      Create file                          */
/*  Chensong Zhang      Oct/15/2021      Format file                          */
/*  Chensong Zhang      Jan/09/2022      Update Doxygen                       */
/*----------------------------------------------------------------------------*/