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

clean up code

Browse Source
This commit is contained in:
James E McClure
2018-10-23 00:50:07 -04:00
parent e6e3316564
commit 2b23951bce
1 changed files with 25 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

78
models/ColorModel.cpp
View File

@@ -539,7 +539,6 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){
// Basic algorithm to
// 1. Copy phase field to CPU
//if (rank==0) printf("MorphInit: copy data \n");
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
double count,count_global,volume_initial,volume_final;
@@ -547,26 +546,17 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm);
volume_initial = count_global;
//if (rank==0) printf("MorphInit: get blob ids \n");
// 2. Identify connected components of phase field -> phase_label
BlobIDstruct new_index;
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
MPI_Barrier(comm);
/* FILE *IDFILE;
sprintf(LocalRankFilename,"Label.%05i.raw",rank);
IDFILE = fopen(LocalRankFilename,"wb");
fwrite(phase_label.data(),4,N,IDFILE);
fclose(IDFILE);
*/
//if (rank==0) printf("MorphInit: label largest feature \n");
// only operate on component "0"
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
@@ -577,43 +567,40 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){
}
}
}
//if (rank==0) printf("MorphInit: generate distance map \n");
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*Dm);
double temp,value;
double factor=0.5/beta;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 3.f ){
value = phase(i,j,k);
if (value > 1.f) value=1.f;
if (value < -1.f) value=-1.f;
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
temp = -factor*log((1.0+value)/(1.0-value));
/// use this approximation close to the object
if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
if (phase_distance(i,j,k) < 3.f ){
value = phase(i,j,k);
if (value > 1.f) value=1.f;
if (value < -1.f) value=-1.f;
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
temp = -factor*log((1.0+value)/(1.0-value));
/// use this approximation close to the object
if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
}
}
}
}
}
// 4. Apply erosion / dilation operation to phase_distance
//if (rank==0) printf("MorphInit: morphological operation \n");
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double walldist=Averages->SDs(i,j,k);
double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f)));
phase_distance(i,j,k) -= wallweight*morph_delta;
}
}
}
//if (rank==0) printf("MorphInit: reinitialize phase field \n");
// 4. Apply erosion / dilation operation to phase_distance
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double walldist=Averages->SDs(i,j,k);
double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f)));
phase_distance(i,j,k) -= wallweight*morph_delta;
}
}
}
// 5. Update phase indicator field based on new distnace
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
@@ -634,35 +621,20 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm);
volume_final=count_global;
if (rank == 0) printf("Morphological operation change volume fraction by %f \n", (volume_final-volume_initial)/volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", (volume_final-volume_initial)/volume_initial);
/* FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(phase.data(),8,N,OUTFILE);
fclose(OUTFILE);
FILE *DISTFILE;
sprintf(LocalRankFilename,"Distance.%05i.raw",rank);
DISTFILE = fopen(LocalRankFilename,"wb");
fwrite(phase_distance.data(),8,N,DISTFILE);
fclose(DISTFILE);
*/
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
// 7. Re-initialize phase field and density
if (rank==0) printf("MorphInit: re-initialize LBM \n");
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){