OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

write id in vis

Browse Source
This commit is contained in:
JamesEMcclure 2019-07-12 09:21:16 -04:00
parent 62a9468c93
commit 7442a31784
4 changed files with 140 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

95
analysis/SubPhase.cpp
View File

@ -624,4 +624,99 @@ void SubPhase::Full(){
}
void SubPhase::AggregateLabels(char *FILENAME){
int nx = Dm->Nx;
int ny = Dm->Ny;
int nz = Dm->Nz;
int npx = Dm->nprocx();
int npy = Dm->nprocy();
int npz = Dm->nprocz();
int ipx = Dm->iproc();
int ipy = Dm->jproc();
int ipz = Dm->kproc();
int full_nx = npx*(nx-2);
int full_ny = npy*(ny-2);
int full_nz = npz*(nz-2);
int local_size = (nx-2)*(ny-2)*(nz-2);
long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
signed char *LocalID;
LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int n = k*nx*ny+j*nx+i;
signed char local_id_val = Dm->id[n];
if (local_id_val > 0){
double value = Phi(i,j,k);
if (value > 0.0) local_id_val = 1;
else local_id_val = 2;
}
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}
}
}
MPI_Barrier(Dm->Comm);
// populate the FullID
if (Dm->rank() == 0){
signed char *FullID;
FullID = new signed char [full_size];
// first handle local ID for rank 0
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int x = i-1;
int y = j-1;
int z = k-1;
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
// next get the local ID from the other ranks
for (int rnk = 1; rnk<Dm->rank(); rnk++){
ipz = rnk / (npx*npy);
ipy = (rnk - ipz*npx*npy) / npx;
ipx = (rnk - ipz*npx*npy - ipy*npx);
int tag = 15+rnk;
MPI_Recv(FullID,local_size,MPI_CHAR,rnk,tag,Dm->Comm,MPI_STATUS_IGNORE);
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int x = i-1 + ipx*(nx-2);
int y = j-1 + ipy*(ny-2);
int z = k-1 + ipz*(nz-2);
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
}
// write the output
FILE *OUTFILE;
OUTFILE = fopen(FILENAME,"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
}
else{
// send LocalID to rank=0
int tag = 15+ Dm->rank();
int dstrank = 0;
MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Dm->Comm);
}
MPI_Barrier(Dm->Comm);
}

3
analysis/SubPhase.h
View File

@ -101,7 +101,8 @@ public:
void Basic();
void Full();
void Write(int time);
void AggregateLabels(char *FILENAME);
private:
FILE *TIMELOG;
FILE *SUBPHASE;

4
analysis/runAnalysis.cpp
View File

@ -232,6 +232,9 @@ public:
fillData.copy(Averages.Vel_z,VelzData);
fillData.copy(Averages.morph_n->label,BlobData);
IO::writeData( timestep, visData, comm.comm );
char CurrentIDFilename[40];
sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
Averages.AggregateLabels(CurrentIDFilename);
PROFILE_STOP("Save Vis",1);
};
private:
@ -1005,4 +1008,3 @@ void runAnalysis::WriteVisData( int timestep, SubPhase &Averages, const double *
PROFILE_STOP("write vis");
}

143
models/ColorModel.cpp
View File

@ -633,6 +633,7 @@ void ScaLBL_ColorModel::Run(){
double volA = Averages->gnb.V;
volA /= Dm->Volume;
volB /= Dm->Volume;;
initial_volume = volA*Dm->Volume;
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
@ -666,10 +667,9 @@ void ScaLBL_ColorModel::Run(){
isSteady = true;
if ( isSteady ){
initial_volume = volA*Dm->Volume;
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
delta_volume_target = Dm->Volume*morph_delta; //*volA ???? // set target volume change
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
Averages->Full();
Averages->Write(timestep);
analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
@ -748,8 +748,7 @@ void ScaLBL_ColorModel::Run(){
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
double effective_seed_water = seed_water*(1.0+volB/volA);
double massChange = SeedPhaseField(effective_seed_water);
double massChange = SeedPhaseField(seed_water);
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
}
else if (USE_MORPHOPEN_OIL){
@ -968,16 +967,13 @@ 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, *Vel_tmp;
double SEED_THRESHOLD = -0.9;
double *Aq_tmp, *Bq_tmp;
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, Velocity, 3*Np*sizeof(double));
/* for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
@ -999,113 +995,54 @@ 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_tmp[n];
uy = Vel_tmp[n+Np];
uz = Vel_tmp[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_tmp[n];
uy = Vel_tmp[n+Np];
uz = Vel_tmp[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++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
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;
// bias toward seed higher flow rate voxels
ux = Vel_tmp[n];
uy = Vel_tmp[n+Np];
uz = Vel_tmp[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;
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] = 0.3333333333333333*water_value;
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
mass_loss += oil_value - dA;
count++;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
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 random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
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;
// bias toward seed higher flow rate voxels
ux = Vel_tmp[n];
uy = Vel_tmp[n+Np];
uz = Vel_tmp[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;
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] = 0.3333333333333333*water_value;
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
mass_loss += oil_value - dA;
count++;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
count= sumReduce( Dm->Comm, count);