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

add image sequence to FlowAdaptor

Browse Source
This commit is contained in:
JamesEMcclure 2021-06-14 16:33:24 -04:00
parent fddd7ae261
commit 7072ede5a4
3 changed files with 336 additions and 301 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
models/ColorModel.cpp
View File

@ -1990,6 +1990,53 @@ FlowAdaptor::~FlowAdaptor(){
}
double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename){
int rank = M.rank;
int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz;
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
M.Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) M.id[i] = M.Mask->id[i]; // save what was read
for (int i=0; i<Nx*Ny*Nz; i++) M.Dm->id[i] = M.Mask->id[i]; // save what was read
double *PhaseLabel;
PhaseLabel = new double[Nx*Ny*Nz];
M.AssignComponentLabels(PhaseLabel);
double Count = 0.0;
double PoreCount = 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++){
if (M.id[Nx*Ny*k+Nx*j+i] == 2){
PoreCount++;
Count++;
}
else if (M.id[Nx*Ny*k+Nx*j+i] == 1){
PoreCount++;
}
}
}
}
Count=M.Dm->Comm.sumReduce( Count);
PoreCount=M.Dm->Comm.sumReduce( PoreCount);
if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_CopyToDevice(M.Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
M.Dm->Comm.barrier();
ScaLBL_D3Q19_Init(M.fq, M.Np);
ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np);
ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np);
M.Dm->Comm.barrier();
ScaLBL_CopyToHost(M.Averages->Phi.data(),M.Phi,Nx*Ny*Nz*sizeof(double));
double saturation = Count/PoreCount;
return saturation;
}
double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
double MASS_FRACTION_CHANGE = 0.0005;
@ -2034,7 +2081,6 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
vax = vay = vaz = 0.0;
vbx = vby = vbz = 0.0;
mass_a = mass_b = 0.0;
double maxSpeed = 0.0;
double localMaxSpeed = 0.0;
for (int k=1; k<Nz-1; k++){
@ -2068,43 +2114,6 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
}
}
maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
/* for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
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];
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];
phi = (dA - dB) / (dA + dB);
Phase[n] = phi;
mass_a += dA;
mass_b += dB;
if (phi > 0.0){
vax += Vel_x[n];
vay += Vel_y[n];
vaz += Vel_z[n];
}
else {
vbx += Vel_x[n];
vby += Vel_y[n];
vbz += Vel_z[n];
}
}
for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
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];
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];
phi = (dA - dB) / (dA + dB);
Phase[n] = phi;
mass_a += dA;
mass_b += dB;
if (phi > 0.0){
vax += Vel_x[n];
vay += Vel_y[n];
vaz += Vel_z[n];
}
else {
vbx += Vel_x[n];
vby += Vel_y[n];
vbz += Vel_z[n];
}
}
*/
mass_a_global = M.Dm->Comm.sumReduce(mass_a);
mass_b_global = M.Dm->Comm.sumReduce(mass_b);
vax_global = M.Dm->Comm.sumReduce(vax);

4
models/ColorModel.h
View File

@ -73,6 +73,8 @@ public:
double *Velocity;
double *Pressure;
void AssignComponentLabels(double *phase);
private:
Utilities::MPI comm;
@ -85,7 +87,6 @@ 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);
@ -97,6 +98,7 @@ public:
FlowAdaptor(ScaLBL_ColorModel &M);
~FlowAdaptor();
double MoveInterface(ScaLBL_ColorModel &M);
double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;

32
tests/lbpm_color_simulator.cpp
View File

@ -66,12 +66,19 @@ int main( int argc, char **argv )
// structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
if (SimulationMode == "development"){
double MLUPS=0.0;
int timestep = 0;
bool ContinueSimulation = true;
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
/* flow adaptor keys to control */
/* Variables for simulation protocols */
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "none" );
/* image sequence protocol */
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
/* flow adaptor keys to control behavior */
int SKIP_TIMESTEPS = 0;
int MAX_STEADY_TIME = 1000000;
double ENDPOINT_THRESHOLD = 0.1;
@ -83,6 +90,8 @@ int main( int argc, char **argv )
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
}
/* analysis keys*/
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
}
@ -97,14 +106,13 @@ int main( int argc, char **argv )
if (ColorModel.timestep > ColorModel.timestepMax){
ContinueSimulation = false;
}
/* update the fractional flow by adding mass */
/* update the fluid configuration with the flow adapter */
int skip_time = 0;
timestep = ColorModel.timestep;
double SaturationChange = 0.0;
double volB = ColorModel.Averages->gwb.V;
double volA = ColorModel.Averages->gnb.V;
double initialSaturation = volB/(volA + volB);
double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M;
double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M;
double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M;
@ -117,7 +125,23 @@ int main( int argc, char **argv )
if (ContinueSimulation){
while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
timestep += ANALYSIS_INTERVAL;
if (PROTOCOL == "fractional flow") {
Adapt.UpdateFractionalFlow(ColorModel);
}
else if (PROTOCOL == "image sequence"){
// Use image sequence
IMAGE_INDEX++;
if (IMAGE_INDEX < IMAGE_COUNT){
std::string next_image = ImageList[IMAGE_INDEX];
if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
ColorModel.color_db->putScalar<int>("image_index",IMAGE_INDEX);
Adapt.ImageInit(ColorModel, next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
ColorModel.timestep = ColorModel.timestepMax;
}
}
MLUPS = ColorModel.Run(timestep);
double volB = ColorModel.Averages->gwb.V;
double volA = ColorModel.Averages->gnb.V;