greyscale model with simple analysis tool
This commit is contained in:
JamesEMcclure
parent
896a3617b8
commit
390556ceb0
@ -1,13 +1,15 @@
|
||||
#include "analysis/SubPhase.h"
|
||||
#include "analysis/GreyPhase.h"
|
||||
|
||||
// Constructor
|
||||
GreyPhase::GreyPhase(std::shared_ptr <Domain> dm):
|
||||
GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr <Domain> dm):
|
||||
Dm(dm)
|
||||
{
|
||||
Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
|
||||
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0;
|
||||
|
||||
// Global arrays
|
||||
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
|
||||
Porosity.resize(Nx,Ny,Nz); Porosity.fill(0);
|
||||
PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0);
|
||||
Rho_n.resize(Nx,Ny,Nz); Rho_n.fill(0);
|
||||
Rho_w.resize(Nx,Ny,Nz); Rho_w.fill(0);
|
||||
@ -39,17 +41,17 @@ GreyPhase::GreyPhase(std::shared_ptr <Domain> dm):
|
||||
|
||||
|
||||
// Destructor
|
||||
GreyPhase::~GreyPhase()
|
||||
GreyPhaseAnalysis::~GreyPhaseAnalysis()
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
void GreyPhase::Write(int timestep)
|
||||
void GreyPhaseAnalysis::Write(int timestep)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
void GreyPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double B)
|
||||
void GreyPhaseAnalysis::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double B)
|
||||
{
|
||||
Fx = force_x;
|
||||
Fy = force_y;
|
||||
@ -62,7 +64,7 @@ void GreyPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, do
|
||||
beta = B;
|
||||
}
|
||||
|
||||
void GreyPhase::Basic(){
|
||||
void GreyPhaseAnalysis::Basic(){
|
||||
int i,j,k,n,imin,jmin,kmin,kmax;
|
||||
|
||||
// If external boundary conditions are set, do not average over the inlet
|
||||
@ -72,22 +74,8 @@ void GreyPhase::Basic(){
|
||||
double count_w = 0.0;
|
||||
double count_n = 0.0;
|
||||
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
n = k*Nx*Ny + j*Nx + i;
|
||||
// Compute volume averages
|
||||
if ( Dm->id[n] > 0 ){
|
||||
// compute density
|
||||
double nA = Rho_n(n);
|
||||
double nB = Rho_w(n);
|
||||
double phi = (nA-nB)/(nA+nB);
|
||||
Phi(n) = phi;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Water_local.reset();
|
||||
Oil_local.reset();
|
||||
for (k=kmin; k<kmax; k++){
|
||||
for (j=jmin; j<Ny-1; j++){
|
||||
for (i=imin; i<Nx-1; i++){
|
||||
@ -95,75 +83,56 @@ void GreyPhase::Basic(){
|
||||
// Compute volume averages
|
||||
if ( Dm->id[n] > 0 ){
|
||||
// compute density
|
||||
double nA = Rho_n(n);
|
||||
double nB = Rho_w(n);
|
||||
double phi = (nA-nB)/(nA+nB);
|
||||
/* if ( phi > 0.0 ){
|
||||
nA = 1.0;
|
||||
nb.V += 1.0;
|
||||
nb.M += nA*rho_n;
|
||||
// velocity
|
||||
nb.Px += rho_n*nA*Vel_x(n);
|
||||
nb.Py += rho_n*nA*Vel_y(n);
|
||||
nb.Pz += rho_n*nA*Vel_z(n);
|
||||
}
|
||||
else{
|
||||
nB = 1.0;
|
||||
wb.M += nB*rho_w;
|
||||
wb.V += 1.0;
|
||||
double pressure = Pressure(n);
|
||||
double rho_n = Rho_n(n);
|
||||
double rho_w = Rho_w(n);
|
||||
double porosity = Porosity(n);
|
||||
double vel_x = Vel_x(n);
|
||||
double vel_y = Vel_y(n);
|
||||
double vel_z = Vel_z(n);
|
||||
Water_local.p += pressure*(rho_w)/(rho_n+rho_w);
|
||||
Water_local.M += rho_w*porosity;
|
||||
Water_local.Px += rho_w*vel_x;
|
||||
Water_local.Py += rho_w*vel_y;
|
||||
Water_local.Pz += rho_w*vel_z;
|
||||
Oil_local.p += pressure*(rho_n)/(rho_n+rho_w);
|
||||
Oil_local.M += rho_n*porosity;
|
||||
Oil_local.Px += rho_n*vel_x;
|
||||
Oil_local.Py += rho_n*vel_y;
|
||||
Oil_local.Pz += rho_n*vel_z;
|
||||
|
||||
// velocity
|
||||
wb.Px += rho_w*nB*Vel_x(n);
|
||||
wb.Py += rho_w*nB*Vel_y(n);
|
||||
wb.Pz += rho_w*nB*Vel_z(n);
|
||||
}
|
||||
if ( phi > 0.99 ){
|
||||
nb.p += Pressure(n);
|
||||
count_n += 1.0;
|
||||
}
|
||||
else if ( phi < -0.99 ){
|
||||
wb.p += Pressure(n);
|
||||
count_w += 1.0;
|
||||
}
|
||||
*/
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
/* gwb.V=sumReduce( Dm->Comm, wb.V);
|
||||
gnb.V=sumReduce( Dm->Comm, nb.V);
|
||||
gwb.M=sumReduce( Dm->Comm, wb.M);
|
||||
gnb.M=sumReduce( Dm->Comm, nb.M);
|
||||
gwb.Px=sumReduce( Dm->Comm, wb.Px);
|
||||
gwb.Py=sumReduce( Dm->Comm, wb.Py);
|
||||
gwb.Pz=sumReduce( Dm->Comm, wb.Pz);
|
||||
gnb.Px=sumReduce( Dm->Comm, nb.Px);
|
||||
gnb.Py=sumReduce( Dm->Comm, nb.Py);
|
||||
gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
|
||||
|
||||
count_w=sumReduce( Dm->Comm, count_w);
|
||||
count_n=sumReduce( Dm->Comm, count_n);
|
||||
if (count_w > 0.0)
|
||||
gwb.p=sumReduce( Dm->Comm, wb.p) / count_w;
|
||||
else
|
||||
gwb.p = 0.0;
|
||||
if (count_n > 0.0)
|
||||
gnb.p=sumReduce( Dm->Comm, nb.p) / count_n;
|
||||
else
|
||||
gnb.p = 0.0;
|
||||
Oil.M=sumReduce( Dm->Comm, Oil_local.M);
|
||||
Oil.p=sumReduce( Dm->Comm, Oil_local.p);
|
||||
Oil.Px=sumReduce( Dm->Comm, Oil_local.Px);
|
||||
Oil.Py=sumReduce( Dm->Comm, Oil_local.Py);
|
||||
Oil.Pz=sumReduce( Dm->Comm, Oil_local.Pz);
|
||||
|
||||
Water.M=sumReduce( Dm->Comm, Water_local.M);
|
||||
Water.p=sumReduce( Dm->Comm, Water_local.p);
|
||||
Water.Px=sumReduce( Dm->Comm, Water_local.Px);
|
||||
Water.Py=sumReduce( Dm->Comm, Water_local.Py);
|
||||
Water.Pz=sumReduce( Dm->Comm, Water_local.Pz);
|
||||
|
||||
Oil.p /= Oil.M;
|
||||
Water.p /= Water.M;
|
||||
|
||||
// check for NaN
|
||||
bool err=false;
|
||||
if (gwb.V != gwb.V) err=true;
|
||||
if (gnb.V != gnb.V) err=true;
|
||||
if (gwb.p != gwb.p) err=true;
|
||||
if (gnb.p != gnb.p) err=true;
|
||||
if (gwb.Px != gwb.Px) err=true;
|
||||
if (gwb.Py != gwb.Py) err=true;
|
||||
if (gwb.Pz != gwb.Pz) err=true;
|
||||
if (gnb.Px != gnb.Px) err=true;
|
||||
if (gnb.Py != gnb.Py) err=true;
|
||||
if (gnb.Pz != gnb.Pz) err=true;
|
||||
if (Water.M != Water.M) err=true;
|
||||
if (Water.p != Water.p) err=true;
|
||||
if (Water.Px != Water.Px) err=true;
|
||||
if (Water.Py != Water.Py) err=true;
|
||||
if (Water.Pz != Water.Pz) err=true;
|
||||
|
||||
if (Oil.M != Oil.M) err=true;
|
||||
if (Oil.p != Oil.p) err=true;
|
||||
if (Oil.Px != Oil.Px) err=true;
|
||||
if (Oil.Py != Oil.Py) err=true;
|
||||
if (Oil.Pz != Oil.Pz) err=true;
|
||||
|
||||
if (Dm->rank() == 0){
|
||||
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
|
||||
@ -194,23 +163,21 @@ void GreyPhase::Basic(){
|
||||
dir_z = 1.0;
|
||||
force_mag = 1.0;
|
||||
}
|
||||
double saturation=gwb.V/(gwb.V + gnb.V);
|
||||
double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
|
||||
double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
|
||||
//double total_flow_rate = water_flow_rate + not_water_flow_rate;
|
||||
//double fractional_flow = water_flow_rate / total_flow_rate;
|
||||
saturation=Water.M/(Water.M + Oil.M); // assume constant density
|
||||
water_flow_rate=saturation*(Water.Px*dir_x + Water.Py*dir_y + Water.Pz*dir_z)/Water.M;
|
||||
oil_flow_rate=(1.0-saturation)*(Oil.Px*dir_x + Oil.Py*dir_y + Oil.Pz*dir_z)/Oil.M;
|
||||
|
||||
double h = Dm->voxel_length;
|
||||
double krn = h*h*nu_n*not_water_flow_rate / force_mag ;
|
||||
double krn = h*h*nu_n*oil_flow_rate / force_mag ;
|
||||
double krw = h*h*nu_w*water_flow_rate / force_mag;
|
||||
//printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, gwb.p, gnb.p);
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*water_flow_rate,h*oil_flow_rate, Water.p, Oil.p);
|
||||
fflush(TIMELOG);
|
||||
}
|
||||
*/
|
||||
|
||||
if (err==true){
|
||||
// exception if simulation produceds NaN
|
||||
printf("GreyPhase.cpp: NaN encountered, may need to check simulation parameters \n");
|
||||
printf("GreyPhaseAnalysis.cpp: NaN encountered, may need to check simulation parameters \n");
|
||||
}
|
||||
ASSERT(err==false);
|
||||
}
|
||||
|
||||
@ -6,7 +6,7 @@
|
||||
#define GreyPhase_INC
|
||||
|
||||
#include <vector>
|
||||
#include "common/Domain.h"
|
||||
#include "common/ScaLBL.h"
|
||||
#include "common/Communication.h"
|
||||
#include "analysis/analysis.h"
|
||||
#include "common/Utilities.h"
|
||||
@ -16,6 +16,17 @@
|
||||
#include "IO/Writer.h"
|
||||
|
||||
class GreyPhase{
|
||||
public:
|
||||
double p;
|
||||
double M,Px,Py,Pz;
|
||||
void reset(){
|
||||
p=M=Px=Py=Pz=0.0;
|
||||
}
|
||||
|
||||
private:
|
||||
};
|
||||
|
||||
class GreyPhaseAnalysis{
|
||||
public:
|
||||
std::shared_ptr <Domain> Dm;
|
||||
double Volume;
|
||||
@ -24,9 +35,16 @@ public:
|
||||
double nu_n, nu_w;
|
||||
double gamma_wn, beta;
|
||||
double Fx, Fy, Fz;
|
||||
// outputs
|
||||
double saturation,water_flow_rate, oil_flow_rate;
|
||||
|
||||
//simulation outputs (averaged values)
|
||||
GreyPhase Water, Oil;
|
||||
GreyPhase Water_local, Oil_local;
|
||||
//...........................................................................
|
||||
int Nx,Ny,Nz;
|
||||
int Nx,Ny,Nz;
|
||||
IntArray SDs; // contains porosity map
|
||||
IntArray Porosity; // contains porosity map
|
||||
IntArray PhaseID; // Phase ID array
|
||||
DoubleArray Rho_n; // density field
|
||||
DoubleArray Rho_w; // density field
|
||||
@ -37,8 +55,8 @@ public:
|
||||
DoubleArray Vel_y;
|
||||
DoubleArray Vel_z;
|
||||
|
||||
GreyPhase(std::shared_ptr <Domain> Dm);
|
||||
~GreyPhase();
|
||||
GreyPhaseAnalysis(std::shared_ptr <Domain> Dm);
|
||||
~GreyPhaseAnalysis();
|
||||
|
||||
void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double beta);
|
||||
void Basic();
|
||||
|
||||
@ -135,7 +135,7 @@ void ScaLBL_GreyscaleColorModel::SetDomain(){
|
||||
N = Nx*Ny*Nz;
|
||||
id = new signed char [N];
|
||||
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
|
||||
Averages = std::shared_ptr<GreyPhase> ( new SubPhase(Dm) ); // TwoPhase analysis object
|
||||
Averages = std::shared_ptr<GreyPhaseAnalysis> ( new GreyPhaseAnalysis(Dm) ); // TwoPhase analysis object
|
||||
MPI_Barrier(comm);
|
||||
Dm->CommInit();
|
||||
MPI_Barrier(comm);
|
||||
@ -831,15 +831,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
USE_SEED = true;
|
||||
USE_MORPH = true;
|
||||
}
|
||||
else if (protocol == "open connected oil"){
|
||||
morph_delta = -0.05;
|
||||
USE_MORPH = true;
|
||||
USE_MORPHOPEN_OIL = true;
|
||||
}
|
||||
else if (protocol == "shell aggregation"){
|
||||
morph_delta = -0.05;
|
||||
USE_MORPH = true;
|
||||
}
|
||||
|
||||
if (greyscaleColor_db->keyExists( "capillary_number" )){
|
||||
capillary_number = greyscaleColor_db->getScalar<double>( "capillary_number" );
|
||||
SET_CAPILLARY_NUMBER=true;
|
||||
@ -868,11 +860,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
morph_interval = analysis_db->getScalar<int>( "morph_interval" );
|
||||
USE_MORPH = true;
|
||||
}
|
||||
if (analysis_db->keyExists( "use_morphopen_oil" )){
|
||||
USE_MORPHOPEN_OIL = analysis_db->getScalar<bool>( "use_morphopen_oil" );
|
||||
if (rank == 0 && USE_MORPHOPEN_OIL) printf("Volume change by morphological opening \n");
|
||||
USE_MORPH = true;
|
||||
}
|
||||
if (analysis_db->keyExists( "tolerance" )){
|
||||
tolerance = analysis_db->getScalar<double>( "tolerance" );
|
||||
}
|
||||
@ -908,20 +895,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
printf(" morph_delta = %f \n",morph_delta);
|
||||
printf(" seed_water = %f \n",seed_water);
|
||||
}
|
||||
else if (protocol == "open connected oil"){
|
||||
printf(" using protocol = open connected oil \n");
|
||||
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
|
||||
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
|
||||
printf(" tolerance = %f \n",tolerance);
|
||||
printf(" morph_delta = %f \n",morph_delta);
|
||||
}
|
||||
else if (protocol == "shell aggregation"){
|
||||
printf(" using protocol = shell aggregation \n");
|
||||
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
|
||||
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
|
||||
printf(" tolerance = %f \n",tolerance);
|
||||
printf(" morph_delta = %f \n",morph_delta);
|
||||
}
|
||||
printf("No. of timesteps: %i \n", timestepMax);
|
||||
fflush(stdout);
|
||||
}
|
||||
@ -938,7 +911,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
//std::shared_ptr<Database> analysis_db;
|
||||
bool Regular = false;
|
||||
auto current_db = db->cloneDatabase();
|
||||
runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
|
||||
//runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
|
||||
//analysis.createThreads( analysis_method, 4 );
|
||||
while (timestep < timestepMax ) {
|
||||
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
|
||||
@ -1076,26 +1049,20 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
|
||||
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){
|
||||
printf("%i %f \n",timestep,din);
|
||||
ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure);
|
||||
ScaLBL_Comm->RegularLayout(Map,&Den[0],Averages->Rho_n);
|
||||
ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w);
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x);
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y);
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z);
|
||||
Averages->Basic();
|
||||
}
|
||||
// Run the analysis
|
||||
analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
|
||||
|
||||
// allow initial ramp-up to get closer to steady state
|
||||
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
|
||||
analysis.finish();
|
||||
//analysis.finish();
|
||||
CURRENT_STEADY_TIMESTEPS += analysis_interval;
|
||||
|
||||
double volB = Averages->gwb.V;
|
||||
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;
|
||||
double vB_x = Averages->gwb.Px/Averages->gwb.M;
|
||||
double vB_y = Averages->gwb.Py/Averages->gwb.M;
|
||||
double vB_z = Averages->gwb.Pz/Averages->gwb.M;
|
||||
double muA = rhoA*(tauA-0.5)/3.f;
|
||||
double muB = rhoB*(tauB-0.5)/3.f;
|
||||
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
|
||||
@ -1109,10 +1076,12 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
dir_z = 1.0;
|
||||
force_mag = 1.0;
|
||||
}
|
||||
double current_saturation = volB/(volA+volB);
|
||||
double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
|
||||
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
|
||||
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
|
||||
double current_saturation = Averages->saturation;
|
||||
double volA = current_saturation*Mask->Porosity()*Mask->Volume;
|
||||
double volB = (1.0-current_saturation)*Mask->Porosity()*Mask->Volume;
|
||||
double flow_rate_A = Averages->oil_flow_rate;
|
||||
double flow_rate_B = Averages->water_flow_rate;
|
||||
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(6.0*alpha);
|
||||
|
||||
if ( morph_timesteps > morph_interval ){
|
||||
|
||||
@ -1165,17 +1134,17 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
volA_prev = volA;
|
||||
//******************************** **/
|
||||
/** compute averages & write data **/
|
||||
Averages->Full();
|
||||
/*Averages->Full();
|
||||
Averages->Write(timestep);
|
||||
analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
|
||||
analysis.finish();
|
||||
|
||||
*/
|
||||
if (rank==0){
|
||||
printf("** WRITE STEADY POINT *** ");
|
||||
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
|
||||
double h = Dm->voxel_length;
|
||||
// pressures
|
||||
double pA = Averages->gnb.p;
|
||||
/*double pA = Averages->gnb.p;
|
||||
double pB = Averages->gwb.p;
|
||||
double pAc = Averages->gnc.p;
|
||||
double pBc = Averages->gwc.p;
|
||||
@ -1231,7 +1200,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
|
||||
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
|
||||
fclose(kr_log_file);
|
||||
|
||||
*/
|
||||
printf(" Measured capillary number %f \n ",Ca);
|
||||
}
|
||||
if (SET_CAPILLARY_NUMBER ){
|
||||
@ -1283,16 +1252,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
double massChange = SeedPhaseField(seed_water);
|
||||
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target);
|
||||
}
|
||||
else if (USE_MORPHOPEN_OIL){
|
||||
delta_volume = volA*Dm->Volume - initial_volume;
|
||||
if (rank==0) printf("***Morphological opening of connected oil, target volume change %f ***\n", delta_volume_target);
|
||||
MorphOpenConnected(delta_volume_target);
|
||||
}
|
||||
else {
|
||||
if (rank==0) printf("***Shell aggregation, target volume change %f ***\n", delta_volume_target);
|
||||
//double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A
|
||||
delta_volume += MorphInit(beta,delta_volume_target-delta_volume);
|
||||
}
|
||||
|
||||
if ( (delta_volume - delta_volume_target)/delta_volume_target > 0.0 ){
|
||||
MORPH_ADAPT = false;
|
||||
@ -1316,7 +1275,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
|
||||
}
|
||||
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
|
||||
}
|
||||
analysis.finish();
|
||||
//analysis.finish();
|
||||
PROFILE_STOP("Loop");
|
||||
PROFILE_SAVE("lbpm_color_simulator",1);
|
||||
//************************************************************************
|
||||
@ -1380,149 +1339,6 @@ double ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
|
||||
return saturation;
|
||||
|
||||
}
|
||||
|
||||
double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_change){
|
||||
|
||||
int nx = Nx;
|
||||
int ny = Ny;
|
||||
int nz = Nz;
|
||||
int n;
|
||||
int N = nx*ny*nz;
|
||||
double volume_change=0.0;
|
||||
|
||||
if (target_volume_change < 0.0){
|
||||
Array<char> id_solid(nx,ny,nz);
|
||||
Array<int> phase_label(nx,ny,nz);
|
||||
DoubleArray distance(Nx,Ny,Nz);
|
||||
DoubleArray phase(nx,ny,nz);
|
||||
signed char *id_connected;
|
||||
id_connected = new signed char [nx*ny*nz];
|
||||
|
||||
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
|
||||
|
||||
// Extract only the connected part of NWP
|
||||
BlobIDstruct new_index;
|
||||
double vF=0.0; double vS=0.0;
|
||||
ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
|
||||
MPI_Barrier(Dm->Comm);
|
||||
|
||||
long long count_connected=0;
|
||||
long long count_porespace=0;
|
||||
long long count_water=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++){
|
||||
n=k*nx*ny+j*nx+i;
|
||||
// only apply opening to connected component
|
||||
if ( phase_label(i,j,k) == 0){
|
||||
count_connected++;
|
||||
}
|
||||
if (id[n] > 0){
|
||||
count_porespace++;
|
||||
}
|
||||
if (id[n] == 2){
|
||||
count_water++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
count_connected=sumReduce( Dm->Comm, count_connected);
|
||||
count_porespace=sumReduce( Dm->Comm, count_porespace);
|
||||
count_water=sumReduce( Dm->Comm, count_water);
|
||||
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
n=k*nx*ny+j*nx+i;
|
||||
// only apply opening to connected component
|
||||
if ( phase_label(i,j,k) == 0){
|
||||
id_solid(i,j,k) = 1;
|
||||
id_connected[n] = 2;
|
||||
id[n] = 2;
|
||||
/* delete the connected component */
|
||||
phase(i,j,k) = -1.0;
|
||||
}
|
||||
else{
|
||||
id_solid(i,j,k) = 0;
|
||||
id_connected[n] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
CalcDist(distance,id_solid,*Dm);
|
||||
|
||||
signed char water=2;
|
||||
signed char notwater=1;
|
||||
double SW=-(target_volume_change)/count_connected;
|
||||
MorphOpen(distance, id_connected, Dm, SW, water, notwater);
|
||||
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
n=k*nx*ny+j*nx+i;
|
||||
// only apply opening to connected component
|
||||
if ( id_connected[n] == 1){
|
||||
id_solid(i,j,k) = 0;
|
||||
}
|
||||
else{
|
||||
id_solid(i,j,k) = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
CalcDist(distance,id_solid,*Dm);
|
||||
|
||||
// re-initialize
|
||||
double beta = 0.95;
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
n=k*nx*ny+j*nx+i;
|
||||
double d = distance(i,j,k);
|
||||
if (Averages->SDs(i,j,k) > 0.f){
|
||||
if (d < 3.f){
|
||||
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int count_morphopen=0.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++){
|
||||
n=k*nx*ny+j*nx+i;
|
||||
// only apply opening to connected component
|
||||
if ( id_connected[n] == 1){
|
||||
count_morphopen++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
count_morphopen=sumReduce( Dm->Comm, count_morphopen);
|
||||
volume_change = double(count_morphopen - count_connected);
|
||||
|
||||
if (rank==0) printf(" opening of connected oil %f \n",volume_change/count_connected);
|
||||
|
||||
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
|
||||
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 == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
|
||||
if (Dm->kproc()==0){
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
|
||||
}
|
||||
if (Dm->kproc() == nprocz-1){
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
return(volume_change);
|
||||
}
|
||||
double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
|
||||
srand(time(NULL));
|
||||
double mass_loss =0.f;
|
||||
@ -1598,219 +1414,6 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
|
||||
return(mass_loss);
|
||||
}
|
||||
|
||||
double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double target_delta_volume){
|
||||
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
|
||||
|
||||
double vF = 0.f;
|
||||
double vS = 0.f;
|
||||
double delta_volume;
|
||||
double WallFactor = 0.0;
|
||||
bool USE_CONNECTED_NWP = false;
|
||||
|
||||
DoubleArray phase(Nx,Ny,Nz);
|
||||
IntArray phase_label(Nx,Ny,Nz);;
|
||||
DoubleArray phase_distance(Nx,Ny,Nz);
|
||||
Array<char> phase_id(Nx,Ny,Nz);
|
||||
fillHalo<double> fillDouble(Dm->Comm,Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1);
|
||||
|
||||
|
||||
// Basic algorithm to
|
||||
// 1. Copy phase field to CPU
|
||||
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
|
||||
|
||||
double count = 0.f;
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
|
||||
}
|
||||
}
|
||||
}
|
||||
double volume_initial = sumReduce( Dm->Comm, count);
|
||||
/*
|
||||
sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
|
||||
FILE *INPUT = fopen(LocalRankFilename,"wb");
|
||||
fwrite(phase.data(),8,N,INPUT);
|
||||
fclose(INPUT);
|
||||
*/
|
||||
// 2. Identify connected components of phase field -> phase_label
|
||||
|
||||
double volume_connected = 0.0;
|
||||
double second_biggest = 0.0;
|
||||
if (USE_CONNECTED_NWP){
|
||||
BlobIDstruct new_index;
|
||||
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
|
||||
MPI_Barrier(Dm->Comm);
|
||||
|
||||
// only operate on component "0"
|
||||
count = 0.0;
|
||||
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
int label = phase_label(i,j,k);
|
||||
if (label == 0 ){
|
||||
phase_id(i,j,k) = 0;
|
||||
count += 1.0;
|
||||
}
|
||||
else
|
||||
phase_id(i,j,k) = 1;
|
||||
if (label == 1 ){
|
||||
second_biggest += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
volume_connected = sumReduce( Dm->Comm, count);
|
||||
second_biggest = sumReduce( Dm->Comm, second_biggest);
|
||||
}
|
||||
else {
|
||||
// use the whole NWP
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
if (Averages->SDs(i,j,k) > 0.f){
|
||||
if (phase(i,j,k) > 0.f ){
|
||||
phase_id(i,j,k) = 0;
|
||||
}
|
||||
else {
|
||||
phase_id(i,j,k) = 1;
|
||||
}
|
||||
}
|
||||
else {
|
||||
phase_id(i,j,k) = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*int reach_x, reach_y, reach_z;
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
}
|
||||
}
|
||||
}*/
|
||||
|
||||
// 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;
|
||||
}
|
||||
// erase the original object
|
||||
phase(i,j,k) = -1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (USE_CONNECTED_NWP){
|
||||
if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
|
||||
// if connected volume is less than 2% just delete the whole thing
|
||||
if (rank==0) printf("Connected region has shrunk! \n");
|
||||
REVERSE_FLOW_DIRECTION = true;
|
||||
}
|
||||
|
||||
/* else{*/
|
||||
if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
|
||||
}
|
||||
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
|
||||
double target_delta_volume_incremental = target_delta_volume;
|
||||
if (fabs(target_delta_volume) > 0.01*volume_initial)
|
||||
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
|
||||
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
|
||||
|
||||
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) < 0.0 ) phase_id(i,j,k) = 0;
|
||||
else phase_id(i,j,k) = 1;
|
||||
//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
|
||||
|
||||
// 5. Update phase indicator field based on new distnace
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
double d = phase_distance(i,j,k);
|
||||
if (Averages->SDs(i,j,k) > 0.f){
|
||||
if (d < 3.f){
|
||||
//phase(i,j,k) = -1.0;
|
||||
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
fillDouble.fill(phase);
|
||||
//}
|
||||
|
||||
count = 0.f;
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f){
|
||||
count+=1.f;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
double volume_final= sumReduce( Dm->Comm, count);
|
||||
|
||||
delta_volume = (volume_final-volume_initial);
|
||||
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
|
||||
if (rank == 0) printf(" new saturation = %f \n", volume_final/(0.238323*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs)));
|
||||
|
||||
// 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));
|
||||
/*
|
||||
sprintf(LocalRankFilename,"dist_final.%05i.raw",rank);
|
||||
FILE *DIST = fopen(LocalRankFilename,"wb");
|
||||
fwrite(phase_distance.data(),8,N,DIST);
|
||||
fclose(DIST);
|
||||
|
||||
sprintf(LocalRankFilename,"phi_final.%05i.raw",rank);
|
||||
FILE *PHI = fopen(LocalRankFilename,"wb");
|
||||
fwrite(phase.data(),8,N,PHI);
|
||||
fclose(PHI);
|
||||
*/
|
||||
// 7. Re-initialize phase field and density
|
||||
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 == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
|
||||
if (Dm->kproc()==0){
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
|
||||
}
|
||||
if (Dm->kproc() == nprocz-1){
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
|
||||
}
|
||||
}
|
||||
return delta_volume;
|
||||
}
|
||||
|
||||
void ScaLBL_GreyscaleColorModel::WriteDebug(){
|
||||
// Copy back final phase indicator field and convert to regular layout
|
||||
DoubleArray PhaseField(Nx,Ny,Nz);
|
||||
|
||||
@ -10,8 +10,7 @@ Implementation of two-fluid greyscale color lattice boltzmann model
|
||||
#include <fstream>
|
||||
|
||||
#include "common/Communication.h"
|
||||
#include "analysis/TwoPhase.h"
|
||||
#include "analysis/runAnalysis.h"
|
||||
#include "analysis/GreyPhase.h"
|
||||
#include "common/MPI_Helpers.h"
|
||||
#include "ProfilerApp.h"
|
||||
#include "threadpool/thread_pool.h"
|
||||
@ -50,7 +49,7 @@ public:
|
||||
std::shared_ptr<Domain> Mask; // this domain is for lbm
|
||||
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
|
||||
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
|
||||
std::shared_ptr<GreyPhase> Averages;
|
||||
std::shared_ptr<GreyPhaseAnalysis> Averages;
|
||||
|
||||
// input database
|
||||
std::shared_ptr<Database> db;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user