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

depcrecate internal color flow adaptor

Browse Source
This commit is contained in:
James McClure
2021-09-19 09:52:06 -04:00
parent 229feb93fb
commit 7575d5c116
3 changed files with 69 additions and 734 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
analysis/dcel.h
View File

@@ -4,9 +4,6 @@
#include <vector>
#include "analysis/pmmc.h"
/*
*/
/**
* \class Vertex
* @brief store vertex for DCEL data structure

734
models/ColorModel.cpp
View File

@@ -921,7 +921,11 @@ double ScaLBL_ColorModel::Run(int returntime){
void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
int analysis_interval = 1000; // number of timesteps in between in situ analysis
if (analysis_db->keyExists( "analysis_interval" )){
analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
}
/*
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
@@ -940,7 +944,6 @@ void ScaLBL_ColorModel::Run(){
int CURRENT_MORPH_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
int CURRENT_STEADY_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
int morph_interval = 100000;
int analysis_interval = 1000; // number of timesteps in between in situ analysis
int morph_timesteps = 0;
double morph_delta = 0.0;
double seed_water = 0.0;
@@ -951,7 +954,6 @@ void ScaLBL_ColorModel::Run(){
double delta_volume = 0.0;
double delta_volume_target = 0.0;
/* history for morphological algoirthm */
double KRA_MORPH_FACTOR=0.5;
double volA_prev = 0.0;
double log_krA_prev = 1.0;
@@ -1003,9 +1005,7 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "tolerance" )){
tolerance = analysis_db->getScalar<double>( "tolerance" );
}
if (analysis_db->keyExists( "analysis_interval" )){
analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
}
if (analysis_db->keyExists( "min_steady_timesteps" )){
MIN_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "min_steady_timesteps" );
}
@@ -1016,7 +1016,6 @@ void ScaLBL_ColorModel::Run(){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
/* defaults for simulation protocols */
auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){
// Get the list of images
@@ -1062,8 +1061,7 @@ void ScaLBL_ColorModel::Run(){
double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB;
}
}
if (rank==0){
} if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
printf(" using protocol = image sequence \n");
@@ -1118,7 +1116,7 @@ void ScaLBL_ColorModel::Run(){
printf("No. of timesteps: %i \n", timestepMax);
fflush(stdout);
}
*/
//************ MAIN ITERATION LOOP ***************************************/
comm.barrier();
PROFILE_START("Loop");
@@ -1218,242 +1216,6 @@ void ScaLBL_ColorModel::Run(){
}
// 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();
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);
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
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);
if ( morph_timesteps > morph_interval ){
bool isSteady = false;
if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_STEADY_TIMESTEPS > MIN_STEADY_TIMESTEPS))
isSteady = true;
if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
isSteady = true;
if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_STEADY_TIMESTEPS > RESCALE_FORCE_AFTER_TIMESTEP){
RESCALE_FORCE = false;
double RESCALE_FORCE_FACTOR = capillary_number / Ca;
if (RESCALE_FORCE_FACTOR > 2.0) RESCALE_FORCE_FACTOR = 2.0;
if (RESCALE_FORCE_FACTOR < 0.5) RESCALE_FORCE_FACTOR = 0.5;
Fx *= RESCALE_FORCE_FACTOR;
Fy *= RESCALE_FORCE_FACTOR;
Fz *= RESCALE_FORCE_FACTOR;
force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
color_db->putVector<double>("F",{Fx,Fy,Fz});
}
if ( isSteady ){
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
//****** ENDPOINT ADAPTATION ********/
double krA_TMP= fabs(muA*flow_rate_A / force_mag);
double krB_TMP= fabs(muB*flow_rate_B / force_mag);
log_krA = log(krA_TMP);
if (krA_TMP < 0.0){
// cannot do endpoint adaptation if kr is negative
log_krA = log_krA_prev;
}
else if (krA_TMP < krB_TMP && morph_delta > 0.0){
/** morphological target based on relative permeability for A **/
log_krA_target = log(KRA_MORPH_FACTOR*(krA_TMP));
slope_krA_volume = (log_krA - log_krA_prev)/(Dm->Volume*(volA - volA_prev));
delta_volume_target=min(delta_volume_target,Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume));
if (rank==0){
printf(" Enabling endpoint adaptation: krA = %f, krB = %f \n",krA_TMP,krB_TMP);
printf(" log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));
}
}
log_krA_prev = log_krA;
volA_prev = volA;
//******************************** **/
/** compute averages & write data **/
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 pB = Averages->gwb.p;
double pAc = Averages->gnc.p;
double pBc = Averages->gwc.p;
double pAB = (pA-pB)/(h*6.0*alpha);
double pAB_connected = (pAc-pBc)/(h*6.0*alpha);
// connected contribution
double Vol_nc = Averages->gnc.V/Dm->Volume;
double Vol_wc = Averages->gwc.V/Dm->Volume;
double Vol_nd = Averages->gnd.V/Dm->Volume;
double Vol_wd = Averages->gwd.V/Dm->Volume;
double Mass_n = Averages->gnc.M + Averages->gnd.M;
double Mass_w = Averages->gwc.M + Averages->gwd.M;
double vAc_x = Averages->gnc.Px/Mass_n;
double vAc_y = Averages->gnc.Py/Mass_n;
double vAc_z = Averages->gnc.Pz/Mass_n;
double vBc_x = Averages->gwc.Px/Mass_w;
double vBc_y = Averages->gwc.Py/Mass_w;
double vBc_z = Averages->gwc.Pz/Mass_w;
// disconnected contribution
double vAd_x = Averages->gnd.Px/Mass_n;
double vAd_y = Averages->gnd.Py/Mass_n;
double vAd_z = Averages->gnd.Pz/Mass_n;
double vBd_x = Averages->gwd.Px/Mass_w;
double vBd_y = Averages->gwd.Py/Mass_w;
double vBd_z = Averages->gwd.Pz/Mass_w;
double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z);
double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z);
double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z);
double flow_rate_B_disconnected = (Vol_wd)*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z);
double kAeff_connected = h*h*muA*flow_rate_A_connected/(force_mag);
double kBeff_connected = h*h*muB*flow_rate_B_connected/(force_mag);
double kAeff_disconnected = h*h*muA*flow_rate_A_disconnected/(force_mag);
double kBeff_disconnected = h*h*muB*flow_rate_B_disconnected/(force_mag);
double kAeff = h*h*muA*(flow_rate_A)/(force_mag);
double kBeff = h*h*muB*(flow_rate_B)/(force_mag);
double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag;
double Mobility = muA/muB;
bool WriteHeader=false;
FILE * kr_log_file = fopen("relperm.csv","r");
if (kr_log_file != NULL)
fclose(kr_log_file);
else
WriteHeader=true;
kr_log_file = fopen("relperm.csv","a");
if (WriteHeader)
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
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 ){
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
color_db->putVector<double>("F",{Fx,Fy,Fz});
}
CURRENT_STEADY_TIMESTEPS = 0;
}
else{
if (rank==0){
printf("** Continue to simulate steady *** \n ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
}
}
morph_timesteps=0;
Ca_previous = Ca;
}
if (MORPH_ADAPT ){
CURRENT_MORPH_TIMESTEPS += analysis_interval;
if (USE_DIRECT){
// Use image sequence
IMAGE_INDEX++;
MORPH_ADAPT = false;
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);
color_db->putScalar<int>("image_index",IMAGE_INDEX);
ImageInit(next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
timestep = timestepMax;
}
}
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
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;
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
else if (!(USE_DIRECT) && CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
initial_volume = volA*Dm->Volume;
delta_volume = 0.0;
RESCALE_FORCE = true;
if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
RESCALE_FORCE = true;
}
}
morph_timesteps += analysis_interval;
}
comm.barrier();
}
analysis.finish();
PROFILE_STOP("Loop");
@@ -1477,486 +1239,6 @@ void ScaLBL_ColorModel::Run(){
// ************************************************************************
}
double ScaLBL_ColorModel::ImageInit(std::string Filename){
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
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; // save what was read
double *PhaseLabel;
PhaseLabel = new double[Nx*Ny*Nz];
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 (id[Nx*Ny*k+Nx*j+i] == 2){
PoreCount++;
Count++;
}
else if (id[Nx*Ny*k+Nx*j+i] == 1){
PoreCount++;
}
}
}
}
Count=Dm->Comm.sumReduce( Count);
PoreCount=Dm->Comm.sumReduce( PoreCount);
if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
comm.barrier();
ScaLBL_D3Q19_Init(fq, Np);
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);
comm.barrier();
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
double saturation = Count/PoreCount;
return saturation;
}
double ScaLBL_ColorModel::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
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);
comm.barrier();
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=Dm->Comm.sumReduce( count_connected);
count_porespace=Dm->Comm.sumReduce( count_porespace);
count_water=Dm->Comm.sumReduce( 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=Dm->Comm.sumReduce( 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_ColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
double count =0.f;
double *Aq_tmp, *Bq_tmp;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double 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];
double 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];
double phase_id = (dA - dB) / (dA + dB);
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double 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];
double 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];
double phase_id = (dA - dB) / (dA + dB);
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
count= Dm->Comm.sumReduce( count);
mass_loss= Dm->Comm.sumReduce( mass_loss);
if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
ScaLBL_CopyToDevice(Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(Bq, Bq_tmp, 7*Np*sizeof(double));
return(mass_loss);
}
double ScaLBL_ColorModel::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 = 1.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 = Dm->Comm.sumReduce( count);
double PoreVolume = Dm->Volume*Dm->Porosity();
/*ensure target isn't an absurdly small fraction of pore volume */
if (volume_initial < target_delta_volume*PoreVolume){
volume_initial = target_delta_volume*PoreVolume;
}
/*
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){
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
comm.barrier();
// 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 = Dm->Comm.sumReduce( count);
second_biggest = Dm->Comm.sumReduce( 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= Dm->Comm.sumReduce( 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/(Mask->Porosity()*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_ColorModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);

66
models/ColorModel.h
View File

@@ -21,21 +21,78 @@ Implementation of color lattice boltzmann model
#ifndef ScaLBL_ColorModel_INC
#define ScaLBL_ColorModel_INC
/**
* \class ScaLBL_ColorModel
*
* @details
* The ScaLBL_ColorModel class contains routines to initialize and run a two-component color lattice Boltzmann model
* Momentum transport equations are described by a D3Q19 scheme
* Mass transport equations are described by D3Q7 scheme
*/
class ScaLBL_ColorModel{
public:
/**
* \brief Constructor
* @param RANK processor rank
* @param NP number of processors
* @param COMM MPI communicator
*/
ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM);
~ScaLBL_ColorModel();
// functions in they should be run
/**
* \brief Read simulation parameters
* @param filename input database file that includes "Color" section
*/
void ReadParams(string filename);
/**
* \brief Read simulation parameters
* @param db0 input database that includes "Color" section
*/
void ReadParams(std::shared_ptr<Database> db0);
/**
* \brief Create domain data structures
*/
void SetDomain();
/**
* \brief Read image data
*/
void ReadInput();
/**
* \brief Create color model data structures
*/
void Create();
/**
* \brief Initialize the simulation
*/
void Initialize();
/**
* \brief Run the simulation
*/
void Run();
/**
* \brief Run the simulation
* @param returntime - timestep at which the routine will return
*/
double Run(int returntime);
/**
* \brief Debugging function to dump simulation state to disk
*/
void WriteDebug();
/**
* \brief Copy the phase field for use by external methods
* @param f - DoubleArray to hold the phase field
*/
void getPhaseField(DoubleArray &f);
bool Restart,pBC;
@@ -74,6 +131,9 @@ public:
double *Velocity;
double *Pressure;
/**
* \brief Assign wetting affinity values
*/
void AssignComponentLabels(double *phase);
private:
@@ -88,10 +148,6 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
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);
};
#endif