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

Merge branch 'morphLBM' of github.com:JamesEMcClure/LBPM-WIA into morphLBM

Browse Source
This commit is contained in:
JamesEMcclure 2019-08-21 10:02:50 -04:00
commit 8210b4c0bd
10 changed files with 200 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
analysis/SubPhase.cpp
View File

@ -112,9 +112,9 @@ void SubPhase::Write(int timestep)
fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
@ -509,10 +509,6 @@ void SubPhase::Full(){
double vol_wc_bulk = 0.0;
double vol_nd_bulk = 0.0;
double vol_wd_bulk = 0.0;
double count_wc = 0.0;
double count_nc = 0.0;
double count_wd = 0.0;
double count_nd = 0.0;
for (k=kmin; k<kmax; k++){
for (j=jmin; j<Ny-1; j++){
for (i=imin; i<Nx-1; i++){
@ -537,28 +533,27 @@ void SubPhase::Full(){
}
else if ( phi > 0.0){
if (morph_n->label(i,j,k) > 0 ){
count_nd += 1.0;
vol_nd_bulk += 1.0;
nd.p += Pressure(n);
}
else{
count_nc += 1.0;
vol_nc_bulk += 1.0;
nc.p += Pressure(n);
}
}
else{
// water region
if (morph_w->label(i,j,k) > 0 ){
count_wd += 1.0;
vol_wd_bulk += 1.0;
wd.p += Pressure(n);
}
else{
count_wc += 1.0;
vol_wc_bulk += 1.0;
wc.p += Pressure(n);
}
}
if ( phi > 0.0){
if (morph_n->label(i,j,k) > 0 ){
vol_nd_bulk += 1.0;
nd.M += nA*rho_n;
nd.Px += nA*rho_n*ux;
nd.Py += nA*rho_n*uy;
@ -566,7 +561,6 @@ void SubPhase::Full(){
nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
}
else{
vol_nc_bulk += 1.0;
nc.M += nA*rho_n;
nc.Px += nA*rho_n*ux;
nc.Py += nA*rho_n*uy;
@ -577,7 +571,6 @@ void SubPhase::Full(){
else{
// water region
if (morph_w->label(i,j,k) > 0 ){
vol_wd_bulk += 1.0;
wd.M += nB*rho_w;
wd.Px += nB*rho_w*ux;
wd.Py += nB*rho_w*uy;
@ -585,7 +578,6 @@ void SubPhase::Full(){
wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
}
else{
vol_wc_bulk += 1.0;
wc.M += nB*rho_w;
wc.Px += nB*rho_w*ux;
wc.Py += nB*rho_w*uy;
@ -597,14 +589,6 @@ void SubPhase::Full(){
}
}
}
count_wc=sumReduce( Dm->Comm, count_wc);
count_nc=sumReduce( Dm->Comm, count_nc);
count_wd=sumReduce( Dm->Comm, count_wd);
count_nd=sumReduce( Dm->Comm, count_nd);
gnd.p=sumReduce( Dm->Comm, nd.p) / count_nd;
gwd.p=sumReduce( Dm->Comm, wd.p) / count_wd;
gnc.p=sumReduce( Dm->Comm, nc.p) / count_nc;
gwc.p=sumReduce( Dm->Comm, wc.p) / count_wc;
gnd.M=sumReduce( Dm->Comm, nd.M);
gnd.Px=sumReduce( Dm->Comm, nd.Px);

15
analysis/runAnalysis.cpp
View File

@ -491,11 +491,9 @@ runAnalysis::commWrapper runAnalysis::getComm( )
/******************************************************************
* Constructor/Destructors *
******************************************************************/
runAnalysis::runAnalysis( std::shared_ptr<Database> db,
const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
int Np, bool Regular, double beta, IntArray Map ):
runAnalysis::runAnalysis(std::shared_ptr<Database> db, const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
int Np, bool Regular, IntArray Map ):
d_Np( Np ),
d_beta( beta ),
d_regular ( Regular),
d_rank_info( rank_info ),
d_Map( Map ),
@ -533,6 +531,8 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
auto restart_file = db->getScalar<std::string>( "restart_file" );
d_restartFile = restart_file + "." + rankString;
d_rank = MPI_WORLD_RANK();
writeIDMap(ID_map_struct(),0,id_map_filename);
// Initialize IO for silo
@ -729,10 +729,12 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
/******************************************************************
* Run the analysis *
******************************************************************/
void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
void runAnalysis::run( std::shared_ptr<Database> db, TwoPhase& Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den)
{
int N = d_N[0]*d_N[1]*d_N[2];
int timestep = db->getWithDefault<int>( "timestep", 0 );
// Check which analysis steps we need to perform
auto type = computeAnalysisType( timestep );
@ -897,11 +899,12 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
/******************************************************************
* Run the analysis *
******************************************************************/
void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
void runAnalysis::basic( std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
{
int N = d_N[0]*d_N[1]*d_N[2];
// Check which analysis steps we need to perform
int timestep = db->getWithDefault<int>( "timestep", 0 );
auto type = computeAnalysisType( timestep );
if ( type == AnalysisType::AnalyzeNone )
return;

6
analysis/runAnalysis.h
View File

@ -40,16 +40,16 @@ public:
//! Constructor
runAnalysis( std::shared_ptr<Database> db, const RankInfoStruct& rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> dm, int Np, bool Regular, double beta, IntArray Map );
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> dm, int Np, bool Regular, IntArray Map );
//! Destructor
~runAnalysis();
//! Run the next analysis
void run( int timestep, TwoPhase &Averages, const double *Phi,
void run( std::shared_ptr<Database> db, TwoPhase &Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den );
void basic( int timestep, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den );
void basic( std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den );
void WriteVisData( int timestep, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den);
//! Finish all active analysis

38
common/Domain.cpp
View File

@ -42,7 +42,7 @@ static inline void fgetl( char * str, int num, FILE * stream )
********************************************************/
Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
double lx, double ly, double lz, int BC):
Nx(0), Ny(0), Nz(0),
database(NULL), Nx(0), Ny(0), Nz(0),
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), voxel_length(1),
Comm(MPI_COMM_WORLD),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
@ -88,7 +88,7 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
initialize( db );
}
Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
Nx(0), Ny(0), Nz(0),
database(db), Nx(0), Ny(0), Nz(0),
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0),
Comm(MPI_COMM_NULL),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
@ -241,7 +241,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!");
}
void Domain::Decomp(std::shared_ptr<Database> domain_db )
void Domain::Decomp(std::string Filename)
{
//.......................................................................
// Reading the domain information file
@ -266,38 +266,38 @@ void Domain::Decomp(std::shared_ptr<Database> domain_db )
checkerSize = 32;
// Read domain parameters
auto Filename = domain_db->getScalar<std::string>( "Filename" );
//auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto SIZE = domain_db->getVector<int>( "N" );
auto nproc = domain_db->getVector<int>( "nproc" );
if (domain_db->keyExists( "offset" )){
auto offset = domain_db->getVector<int>( "offset" );
//auto Filename = database->getScalar<std::string>( "Filename" );
//auto L = database->getVector<double>( "L" );
auto size = database->getVector<int>( "n" );
auto SIZE = database->getVector<int>( "N" );
auto nproc = database->getVector<int>( "nproc" );
if (database->keyExists( "offset" )){
auto offset = database->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
if (domain_db->keyExists( "InletLayers" )){
auto InletCount = domain_db->getVector<int>( "InletLayers" );
if (database->keyExists( "InletLayers" )){
auto InletCount = database->getVector<int>( "InletLayers" );
inlet_layers_x = InletCount[0];
inlet_layers_y = InletCount[1];
inlet_layers_z = InletCount[2];
}
if (domain_db->keyExists( "OutletLayers" )){
auto OutletCount = domain_db->getVector<int>( "OutletLayers" );
if (database->keyExists( "OutletLayers" )){
auto OutletCount = database->getVector<int>( "OutletLayers" );
outlet_layers_x = OutletCount[0];
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
if (domain_db->keyExists( "checkerSize" )){
checkerSize = domain_db->getScalar<int>( "checkerSize" );
if (database->keyExists( "checkerSize" )){
checkerSize = database->getScalar<int>( "checkerSize" );
}
else {
checkerSize = SIZE[0];
}
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
auto ReadType = domain_db->getScalar<std::string>( "ReadType" );
auto ReadValues = database->getVector<int>( "ReadValues" );
auto WriteValues = database->getVector<int>( "WriteValues" );
auto ReadType = database->getScalar<std::string>( "ReadType" );
if (ReadType == "8bit"){
}

3
common/Domain.h
View File

@ -122,6 +122,7 @@ private:
public: // Public variables (need to create accessors instead)
std::shared_ptr<Database> database;
double Lx,Ly,Lz,Volume,voxel_length;
int Nx,Ny,Nz,N;
int inlet_layers_x, inlet_layers_y, inlet_layers_z;
@ -190,7 +191,7 @@ public: // Public variables (need to create accessors instead)
signed char *id;
void ReadIDs();
void Decomp(std::shared_ptr<Database> domain_db);
void Decomp(std::string Filename);
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();

186
models/ColorModel.cpp
View File

@ -167,19 +167,27 @@ void ScaLBL_ColorModel::SetDomain(){
void ScaLBL_ColorModel::ReadInput(){
if (domain_db->keyExists( "Filename" )){
Mask->Decomp(domain_db);
}
else{
size_t readID;
Mask->ReadIDs();
}
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
if (color_db->keyExists( "image_sequence" )){
auto ImageList = color_db->getVector<std::string>( "image_sequence");
int IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
int IMAGE_COUNT = ImageList.size();
std::string first_image = ImageList[IMAGE_INDEX];
Mask->Decomp(first_image);
IMAGE_INDEX++;
}
else if (color_db->keyExists( "Filename" )){
auto Filename = domain_db->getScalar<std::string>( "Filename" );
Mask->Decomp(Filename);
}
else{
Mask->ReadIDs();
}
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz);
@ -468,6 +476,10 @@ void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
bool SET_CAPILLARY_NUMBER = false;
bool MORPH_ADAPT = false;
bool USE_MORPH = false;
@ -493,6 +505,31 @@ void ScaLBL_ColorModel::Run(){
double delta_volume_target = 0.0;
double RESIDUAL_ENDPOINT_THRESHOLD = 0.04;
auto protocol = color_db->getWithDefault<std::string>( "protocol", "default" );
if (protocol == "image sequence"){
// Get the list of images
USE_DIRECT = true;
ImageList = color_db->getVector<std::string>( "image_sequence");
IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
IMAGE_COUNT = ImageList.size();
IMAGE_INDEX++; // first image is already loaded as initial condition
}
else if (protocol == "seed water"){
morph_delta = 0.05;
seed_water = 0.01;
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 (color_db->keyExists( "residual_endpoint_threshold" )){
RESIDUAL_ENDPOINT_THRESHOLD = color_db->getScalar<double>( "residual_endpoint_threshold" );
}
@ -517,10 +554,6 @@ void ScaLBL_ColorModel::Run(){
morph_interval = analysis_db->getScalar<int>( "morph_interval" );
USE_MORPH = true;
}
if (analysis_db->keyExists( "use_direct" )){
USE_DIRECT = analysis_db->getScalar<bool>( "use_direct" );
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");
@ -541,23 +574,56 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "max_morph_timesteps" )){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
printf(" using protocol = image sequence \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
printf(" tolerance = %f \n",tolerance);
std::string first_image = ImageList[IMAGE_INDEX];
printf(" first image in sequence: %s ***\n", first_image.c_str());
}
else if (protocol == "seed water"){
printf(" using protocol = seed water \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(" 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);
}
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier();
MPI_Barrier(comm);
starttime = MPI_Wtime();
//.........................................
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
//std::shared_ptr<Database> analysis_db;
bool Regular = false;
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, beta, Map );
runAnalysis analysis( analysis_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)); }
@ -639,7 +705,9 @@ void ScaLBL_ColorModel::Run(){
// Run the analysis
//analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.basic( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis_db->putScalar<int>("timestep",timestep);
analysis.basic( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
printf("....inlet pressure=%f \n",din);
@ -757,14 +825,17 @@ void ScaLBL_ColorModel::Run(){
if (MORPH_ADAPT ){
CURRENT_MORPH_TIMESTEPS += analysis_interval;
if (USE_DIRECT){
delta_volume = volA*Dm->Volume - initial_volume;
BoundaryCondition = 4;
flux = morph_delta*5.796*alpha*capillary_number/(fabs(morph_delta)*muA*Dm->Lz);
Fx = Fy = Fz = 0.0;
ScaLBL_Comm->BoundaryCondition = 4;
ScaLBL_Comm_Regular->BoundaryCondition = 4;
if (rank==0) printf("***Direct simulation with flux(A) %f, volume change %f / %f ***\n",flux, delta_volume, delta_volume_target);
// Use image sequence
std::string next_image = ImageList[IMAGE_INDEX];
if (IMAGE_INDEX < IMAGE_COUNT){
if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
ImageInit(next_image);
IMAGE_INDEX++;
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
timestep = timestepMax;
}
}
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
@ -774,11 +845,11 @@ void ScaLBL_ColorModel::Run(){
}
else if (USE_MORPHOPEN_OIL){
delta_volume = volA*Dm->Volume - initial_volume;
if (rank==0) printf("***Morphological opening of connected oil, with target volume change %f ***\n", delta_volume_target);
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("***Morphological step with target volume change %f ***\n", delta_volume_target);
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);
}
@ -789,12 +860,12 @@ void ScaLBL_ColorModel::Run(){
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
if (USE_DIRECT){
BoundaryCondition = 0;
ScaLBL_Comm->BoundaryCondition = 0;
ScaLBL_Comm_Regular->BoundaryCondition = 0;
Fx = capillary_number*dir_x*force_mag / Ca;
Fy = capillary_number*dir_y*force_mag / Ca;
Fz = capillary_number*dir_z*force_mag / Ca;
//BoundaryCondition = 0;
//ScaLBL_Comm->BoundaryCondition = 0;
//ScaLBL_Comm_Regular->BoundaryCondition = 0;
//Fx = capillary_number*dir_x*force_mag / Ca;
//Fy = capillary_number*dir_y*force_mag / Ca;
//Fz = capillary_number*dir_z*force_mag / Ca;
}
}
else if (!(USE_DIRECT) && CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
@ -840,6 +911,57 @@ void ScaLBL_ColorModel::Run(){
// ************************************************************************
}
double ScaLBL_ColorModel::ImageInit(std::string Filename){
bool suppress = false;
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
double *PhaseLabel;
PhaseLabel = new double[Nx*Ny*Nz];
AssignComponentLabels(PhaseLabel);
// consistency check
double Count = 0.0;
double PoreCount = 0.0;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double distance = Averages->SDs(i,j,k);
if (distance > 0.0){
if (id[Nx*Ny*k+Nx*j+i] == 2){
PoreCount++;
Count++;
}
else if (id[Nx*Ny*k+Nx*j+i] == 1){
PoreCount++;
}
/* else if (suppress == false){
printf("WARNING (ScaLBLColorModel::ImageInit) image input file sequence may not be labeled correctly (rank=%i) \n",rank);
suppress = true;
}
*/
}
}
}
}
Count=sumReduce( Dm->Comm, Count);
PoreCount=sumReduce( Dm->Comm, PoreCount);
if (rank==0) printf(" new saturation: %f \n", Count / PoreCount);
ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
MPI_Barrier(comm);
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);
MPI_Barrier(comm);
double saturation = Count/PoreCount;
return saturation;
}
double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
int nx = Nx;

1
models/ColorModel.h
View File

@ -94,6 +94,7 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignComponentLabels(double *phase);
double ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);
double MorphOpenConnected(double target_volume_change);

4
models/DFHModel.cpp
View File

@ -492,7 +492,7 @@ void ScaLBL_DFHModel::Run(){
bool Regular = true;
PROFILE_START("Loop");
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, beta, Map );
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
while (timestep < timestepMax ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
@ -573,7 +573,7 @@ void ScaLBL_DFHModel::Run(){
PROFILE_STOP("Update");
// Run the analysis
analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.run( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
}
analysis.finish();
PROFILE_STOP("Loop");

4
tests/TestBubbleDFH.cpp
View File

@ -405,7 +405,7 @@ int main(int argc, char **argv)
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
//std::shared_ptr<Database> analysis_db;
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, pBC, beta, Map );
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, pBC, Map );
//analysis.createThreads( analysis_method, 4 );
while (timestep < timestepMax && err > tol ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
@ -487,7 +487,7 @@ int main(int argc, char **argv)
PROFILE_STOP("Update");
// Run the analysis
analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.run( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
}
analysis.finish();

1
tests/TestDatabase.cpp
View File

@ -37,6 +37,7 @@ int main(int argc, char **argv)
// write a simple database and test that it can be read by LBPM
auto db = std::make_shared<Database>( );
db->putScalar<int>( "BC", BC );
db->putScalar<int>( "BC", BC );
db->putVector<int>( "nproc", { npx, npy, npz } );
db->putVector<int>( "n", { nx, ny, nz } );
db->putVector<std::string>( "Files", List);