diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 4a1537d9..9feef866 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -88,12 +88,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np); diff --git a/cuda/GreyscaleColor.cu b/cuda/GreyscaleColor.cu index d5ab9460..f26b1357 100644 --- a/cuda/GreyscaleColor.cu +++ b/cuda/GreyscaleColor.cu @@ -1450,7 +1450,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ @@ -1478,6 +1478,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int double Fx,Fy,Fz; double Fcpx,Fcpy,Fcpz;//capillary penalty force double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1504,6 +1505,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int porosity = Poros[n]; perm = Perm[n]; W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -2165,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2184,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2204,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2226,7 +2241,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -2249,6 +2264,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis double Fx,Fy,Fz; double Fcpx,Fcpy,Fcpz;//capillary penalty force double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -2276,6 +2292,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis porosity = Poros[n]; perm = Perm[n]; W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -2870,6 +2888,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2886,6 +2908,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2901,6 +2927,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -4500,11 +4530,11 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm, Vel, Pressure, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -4515,11 +4545,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm,Vel,Pressure, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index 73c1f878..9e3a9bd5 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -309,18 +309,26 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt // oil-wet < 0 // neutral = 0 (i.e. no penalty) double *GreySolidW_host = new double [Np]; + double *GreySn_host = new double [Np]; + double *GreySw_host = new double [Np]; size_t NLABELS=0; signed char VALUE=0; double AFFINITY=0.f; + double Sn,Sw;//end-point saturation of greynodes set by users auto LabelList = greyscaleColor_db->getVector( "GreySolidLabels" ); auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); + auto SnList = greyscaleColor_db->getVector( "GreySnList" ); + auto SwList = greyscaleColor_db->getVector( "GreySwList" ); NLABELS=LabelList.size(); if (NLABELS != AffinityList.size()){ ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); } + if (NLABELS != SnList.size() || NLABELS != SwList.size()){ + ERROR("Error: GreySolidLabels, GreySnList, and GreySwList must be the same length! \n"); + } for (int k=0;k0: water-wet || grey-solid affinity<0: oil-wet \n"); } ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double)); + ScaLBL_CopyToDevice(GreySn, GreySn_host, Np*sizeof(double)); + ScaLBL_CopyToDevice(GreySw, GreySw_host, Np*sizeof(double)); ScaLBL_Comm->Barrier(); delete [] GreySolidW_host; + delete [] GreySn_host; + delete [] GreySw_host; } ////----------------------------------------------------------------------------------------------------------// @@ -598,6 +619,8 @@ void ScaLBL_GreyscaleColorModel::Create(){ //ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz); //ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &GreySn, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np); //........................................................................... @@ -921,7 +944,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ // Halo exchange for phase field ScaLBL_Comm_Regular->SendHalo(Phi); //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //Model-1&4 @@ -950,7 +973,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); //Model-1&4 @@ -983,7 +1006,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } ScaLBL_Comm_Regular->SendHalo(Phi); //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //Model-1&4 @@ -1012,7 +1035,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); //Model-1&4 diff --git a/models/GreyscaleColorModel.h b/models/GreyscaleColorModel.h index a9a0bfc9..bff186e1 100644 --- a/models/GreyscaleColorModel.h +++ b/models/GreyscaleColorModel.h @@ -68,6 +68,8 @@ public: //double *GreySolidPhi; //Model 2 & 3 //double *GreySolidGrad;//Model 1 & 4 double *GreySolidW; + double *GreySn; + double *GreySw; //double *ColorGrad; double *Velocity; double *Pressure;