seed water based on local flow rate

This commit is contained in:
James McClure
2019-06-22 13:50:42 -04:00
parent 34837698cb
commit 6f213544b4

View File

@@ -968,14 +968,16 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
double count =0.f;
double *Aq_tmp, *Bq_tmp;
double *Aq_tmp, *Bq_tmp, *Vel_tmp;
double SEED_THRESHOLD = -0.9;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
Vel_tmp = new double [3*Np];
ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Vel_tmp, Vel, 3*Np*sizeof(double));
/* for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
@@ -997,6 +999,48 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
}
}
*/
/*
*
* look at dot product between velocity field and
* for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
* ux = Vel[n];
* uy = Vel[n+Np];
* uz = Vel[n+2*Np];
* magforce = sqrt(Fx*Fx+ Fy*Fy + Fz*Fz);
* dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce;
* if ( dotprod > flow_rate_A + flow_rate_B ){
*
* }
* }
*
* for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
* }
*
*/
double count_voxels = 0.0;
double average_flow_rate=0.0;
double ux,uy,uz,dotprod;
double magforce = sqrt(Fx*Fx+ Fy*Fy + Fz*Fz);
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
ux = Vel[n];
uy = Vel[n+Np];
uz = Vel[n+2*Np];
dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce;
average_flow_rate += dotprod;
count_voxels += 1.0;
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
ux = Vel[n];
uy = Vel[n+Np];
uz = Vel[n+2*Np];
dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce;
average_flow_rate += dotprod;
count_voxels += 1.0;
}
count_voxels=sumReduce( Dm->Comm, count_voxels);
average_flow_rate=sumReduce( Dm->Comm, average_flow_rate);
average_flow_rate /= count_voxels;
double oil_value = 0.0;
double water_value = 1.0;
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
@@ -1004,7 +1048,14 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
double random_value = double(rand())/ RAND_MAX;
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
// bias toward seed higher flow rate voxels
ux = Vel[n];
uy = Vel[n+Np];
uz = Vel[n+2*Np];
dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce;
if (dotprod > average_flow_rate && phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
Aq_tmp[n] = 0.3333333333333333*oil_value;
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;
@@ -1030,7 +1081,13 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
double random_value = double(rand())/ RAND_MAX;
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
// bias toward seed higher flow rate voxels
ux = Vel[n];
uy = Vel[n+Np];
uz = Vel[n+2*Np];
dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce;
if (dotprod > average_flow_rate && phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
Aq_tmp[n] = 0.3333333333333333*oil_value;
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;