Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA
This commit is contained in:
JamesEMcclure
commit
c80b74a754
@ -93,7 +93,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
|
||||
{
|
||||
// If timelog is empty, write a short header to list the averages
|
||||
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
|
||||
fprintf(TIMELOG,"sw krw krn qw qn pw pn\n");
|
||||
fprintf(TIMELOG,"sw krw krn vw vn pw pn\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -109,35 +109,35 @@ SubPhase::~SubPhase()
|
||||
void SubPhase::Write(int timestep)
|
||||
{
|
||||
if (Dm->rank()==0){
|
||||
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 ",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);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
|
||||
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.p, gwd.p, gnc.p, gnd.p);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.V, gwc.A, gwc.H, gwc.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gnc.V, gnc.A, gnc.H, gnc.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",giwn.V, giwn.A, giwn.H, giwn.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
|
||||
fflush(SUBPHASE);
|
||||
}
|
||||
else{
|
||||
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 ",wc.p, wd.p, nc.p, nd.p);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.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 ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
|
||||
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.p, wd.p, nc.p, nd.p);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.V, wc.A, wc.H, wc.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",nc.V, nc.A, nc.H, nc.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",iwn.V, iwn.A, iwn.H, iwn.X);
|
||||
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
|
||||
}
|
||||
|
||||
}
|
||||
@ -160,9 +160,6 @@ void SubPhase::Basic(){
|
||||
|
||||
// If external boundary conditions are set, do not average over the inlet
|
||||
kmin=1; kmax=Nz-1;
|
||||
if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
|
||||
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
|
||||
|
||||
imin=jmin=1;
|
||||
// If inlet/outlet layers exist use these as default
|
||||
if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
|
||||
@ -307,7 +304,7 @@ void SubPhase::Basic(){
|
||||
double krn = h*h*nu_n*not_water_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*h*h*water_flow_rate,h*h*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*not_water_flow_rate, gwb.p, gnb.p);
|
||||
fflush(TIMELOG);
|
||||
}
|
||||
if (err==true){
|
||||
|
||||
@ -344,7 +344,6 @@ void Domain::Decomp(std::string Filename)
|
||||
char *SegData = NULL;
|
||||
|
||||
if (RANK==0){
|
||||
|
||||
printf("Input media: %s\n",Filename.c_str());
|
||||
printf("Relabeling %lu values\n",ReadValues.size());
|
||||
for (int idx=0; idx<ReadValues.size(); idx++){
|
||||
@ -500,6 +499,7 @@ void Domain::Decomp(std::string Filename)
|
||||
else{
|
||||
// solid checkers
|
||||
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -31,7 +31,6 @@
|
||||
#include "common/Communication.h"
|
||||
#include "common/Database.h"
|
||||
|
||||
|
||||
class Domain;
|
||||
template<class TYPE> class PatchData;
|
||||
|
||||
|
||||
@ -70,6 +70,8 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
|
||||
auto n = domain.getVector<int>( "n" );
|
||||
int rank = comm_rank(MPI_COMM_WORLD);
|
||||
auto nproc = domain.getVector<int>( "nproc" );
|
||||
auto ReadValues = domain.getVector<int>( "ReadValues" );
|
||||
auto WriteValues = domain.getVector<int>( "WriteValues" );
|
||||
RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );
|
||||
|
||||
// Determine the largest file number to get
|
||||
@ -81,18 +83,37 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
|
||||
Array<uint8_t> data;
|
||||
RankInfoStruct srcRankInfo( rank, Nfx, Nfy, Nfz );
|
||||
if ( srcRankInfo.ix >= 0 ) {
|
||||
auto filename = domain.getScalar<std::string>( "Filename" );
|
||||
char tmp[100];
|
||||
if ( filename.find( "0x_0y_0z.gbd.gz" ) != std::string::npos ) {
|
||||
sprintf( tmp, "%ix_%iy_%iz.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
|
||||
filename = filename.replace( filename.find( "0x_0y_0z.gbd.gz" ), 15, std::string( tmp ) );
|
||||
} else if ( filename.find( "x0_y0_z0.gbd.gz" ) != std::string::npos ) {
|
||||
sprintf( tmp, "x%i_y%i_z%i.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
|
||||
filename = filename.replace( filename.find( "x0_y0_z0.gbd.gz" ), 15, std::string( tmp ) );
|
||||
} else {
|
||||
ERROR( "Invalid name for first file" );
|
||||
}
|
||||
data = readMicroCT( filename );
|
||||
auto filename = domain.getScalar<std::string>( "Filename" );
|
||||
char tmp[100];
|
||||
if ( filename.find( "0x_0y_0z.gbd.gz" ) != std::string::npos ) {
|
||||
sprintf( tmp, "%ix_%iy_%iz.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
|
||||
filename = filename.replace( filename.find( "0x_0y_0z.gbd.gz" ), 15, std::string( tmp ) );
|
||||
} else if ( filename.find( "x0_y0_z0.gbd.gz" ) != std::string::npos ) {
|
||||
sprintf( tmp, "x%i_y%i_z%i.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
|
||||
filename = filename.replace( filename.find( "x0_y0_z0.gbd.gz" ), 15, std::string( tmp ) );
|
||||
} else {
|
||||
ERROR( "Invalid name for first file" );
|
||||
}
|
||||
data = readMicroCT( filename );
|
||||
|
||||
// Relabel the data
|
||||
for (int k = 0; k<1024; k++){
|
||||
for (int j = 0; j<1024; j++){
|
||||
for (int i = 0; i<1024; i++){
|
||||
//n = k*Nfx*Nfy + j*Nfx + i;
|
||||
//char locval = loc_id[n];
|
||||
char locval = data(i,j,k);
|
||||
for (int idx=0; idx<ReadValues.size(); idx++){
|
||||
signed char oldvalue=ReadValues[idx];
|
||||
signed char newvalue=WriteValues[idx];
|
||||
if (locval == oldvalue){
|
||||
data(i,j,k) = newvalue;
|
||||
idx = ReadValues.size();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Redistribute the data
|
||||
|
||||
@ -19,6 +19,7 @@ color lattice boltzmann model
|
||||
#include "models/ColorModel.h"
|
||||
#include "analysis/distance.h"
|
||||
#include "analysis/morphology.h"
|
||||
#include "common/ReadMicroCT.h"
|
||||
#include <stdlib.h>
|
||||
#include <time.h>
|
||||
|
||||
@ -133,7 +134,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){
|
||||
inletB=0.f;
|
||||
outletA=0.f;
|
||||
outletB=1.f;
|
||||
if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
|
||||
//if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
|
||||
|
||||
BoundaryCondition = 0;
|
||||
if (domain_db->keyExists( "BC" )){
|
||||
@ -204,6 +205,10 @@ void ScaLBL_ColorModel::ReadInput(){
|
||||
Mask->Decomp(first_image);
|
||||
IMAGE_INDEX++;
|
||||
}
|
||||
else if (domain_db->keyExists( "GridFile" )){
|
||||
auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
|
||||
for (int i=0; i<Nx*Ny*Nz; i++) Mask->id[i] = input_id(i);
|
||||
}
|
||||
else if (domain_db->keyExists( "Filename" )){
|
||||
auto Filename = domain_db->getScalar<std::string>( "Filename" );
|
||||
Mask->Decomp(Filename);
|
||||
@ -1490,4 +1495,68 @@ void ScaLBL_ColorModel::WriteDebug(){
|
||||
OUTFILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,OUTFILE);
|
||||
fclose(OUTFILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
|
||||
FILE *AFILE;
|
||||
sprintf(LocalRankFilename,"A.%05i.raw",rank);
|
||||
AFILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,AFILE);
|
||||
fclose(AFILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
|
||||
FILE *BFILE;
|
||||
sprintf(LocalRankFilename,"B.%05i.raw",rank);
|
||||
BFILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,BFILE);
|
||||
fclose(BFILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
|
||||
FILE *PFILE;
|
||||
sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
|
||||
PFILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,PFILE);
|
||||
fclose(PFILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
|
||||
FILE *VELX_FILE;
|
||||
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
|
||||
VELX_FILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,VELX_FILE);
|
||||
fclose(VELX_FILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
|
||||
FILE *VELY_FILE;
|
||||
sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
|
||||
VELY_FILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,VELY_FILE);
|
||||
fclose(VELY_FILE);
|
||||
|
||||
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
|
||||
FILE *VELZ_FILE;
|
||||
sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
|
||||
VELZ_FILE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(PhaseField.data(),8,N,VELZ_FILE);
|
||||
fclose(VELZ_FILE);
|
||||
|
||||
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
|
||||
// FILE *CGX_FILE;
|
||||
// sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
|
||||
// CGX_FILE = fopen(LocalRankFilename,"wb");
|
||||
// fwrite(PhaseField.data(),8,N,CGX_FILE);
|
||||
// fclose(CGX_FILE);
|
||||
//
|
||||
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField);
|
||||
// FILE *CGY_FILE;
|
||||
// sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank);
|
||||
// CGY_FILE = fopen(LocalRankFilename,"wb");
|
||||
// fwrite(PhaseField.data(),8,N,CGY_FILE);
|
||||
// fclose(CGY_FILE);
|
||||
//
|
||||
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField);
|
||||
// FILE *CGZ_FILE;
|
||||
// sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank);
|
||||
// CGZ_FILE = fopen(LocalRankFilename,"wb");
|
||||
// fwrite(PhaseField.data(),8,N,CGZ_FILE);
|
||||
// fclose(CGZ_FILE);
|
||||
|
||||
}
|
||||
|
||||
@ -209,9 +209,18 @@ void ScaLBL_MRTModel::Run(){
|
||||
int SIZE=Np*sizeof(double);
|
||||
|
||||
if (rank==0){
|
||||
FILE * log_file = fopen("Permeability.csv","a");
|
||||
fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n");
|
||||
fclose(log_file);
|
||||
bool WriteHeader=false;
|
||||
FILE *log_file = fopen("Permeability.csv","r");
|
||||
if (log_file != NULL)
|
||||
fclose(log_file);
|
||||
else
|
||||
WriteHeader=true;
|
||||
|
||||
if (WriteHeader){
|
||||
log_file = fopen("Permeability.csv","a+");
|
||||
fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n");
|
||||
fclose(log_file);
|
||||
}
|
||||
}
|
||||
|
||||
//.......create and start timer............
|
||||
|
||||
Loading…
Reference in New Issue
Block a user