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

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

Browse Source
This commit is contained in:
JamesEMcclure
2021-09-14 09:19:41 -04:00
parent 2745e932f5 c3300baf1e
commit f2e4c20d2f
48 changed files with 7026 additions and 875 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@@ -90,7 +90,7 @@ CHECK_ENABLE_FLAG( USE_DOXYGEN 1 )
CHECK_ENABLE_FLAG( USE_LATEX 1 )
FILE( MAKE_DIRECTORY "${${PROJ}_INSTALL_DIR}/doc" )
IF ( USE_DOXYGEN )
SET( DOXYFILE_LATEX YES )
SET( DOXYFILE_LATEX NO )
SET( DOXYFILE_IN "${${PROJ}_SOURCE_DIR}/doxygen/Doxyfile.in" )
SET( DOXY_HEADER_FILE "${${PROJ}_SOURCE_DIR}/doxygen/html/header.html" )
SET( DOXY_FOOTER_FILE "${${PROJ}_SOURCE_DIR}/doxygen/html/footer.html" )

2
IO/Reader.h
View File

@@ -87,7 +87,7 @@ std::shared_ptr<IO::Variable> getVariable( const std::string &path, const std::s
* @brief Reformat the variable to match the mesh
* @details This function modifies the dimensions of the array to match the mesh
* @param[in] mesh The underlying mesh
* @param[in/out] variable The variable name to read
* @param[in,out] var The variable name to read
*/
void reformatVariable( const IO::Mesh &mesh, IO::Variable &var );

1
IO/Xdmf.h
View File

@@ -105,6 +105,7 @@ public:
* @param[in] type The element type
* @param[in] NumElements The number of elements
* @param[in] dofMap The connectivity information (type x NumElements)
* @param[in] NumNodes The number of nodes
* @param[in] x The x coordinates or the xy/xyz coordinates
* @param[in] y The y coordinates (may be null)
* @param[in] z The z coordinates (may be null)

6
IO/netcdf.h
View File

@@ -123,6 +123,8 @@ Array<TYPE> getAtt( int fid, const std::string &att );
* @brief Write the dimensions
* @details This function writes the grid dimensions to netcdf.
* @param fid Handle to the open file
* @param names
* @param dims
*/
std::vector<int> defDim(
int fid, const std::vector<std::string> &names, const std::vector<int> &dims );
@@ -132,6 +134,10 @@ std::vector<int> defDim(
* @brief Write a variable
* @details This function writes a variable to netcdf.
* @param fid Handle to the open file
* @param var Variable to read
* @param dimids
* @param data
* @param rank_info
*/
template<class TYPE>
void write( int fid, const std::string &var, const std::vector<int> &dimids,

4
IO/silo.h
View File

@@ -49,7 +49,7 @@ void close( DBfile *fid );
* @param[in] fid Handle to the open file
* @param[in] name Name of variable
*/
DataType varDataType( DBfile *dbfile, const std::string &name );
DataType varDataType( DBfile *fid, const std::string &name );
/*!
@@ -250,8 +250,6 @@ void writeMultiMesh( DBfile *fid, const std::string &meshname,
* @param[in] varname Mesh name
* @param[in] subVarNames Names of the sub meshes in the form "filename:meshname"
* @param[in] subVarTypes Type of each submesh
* @param[in] ndim Dimension of variable (used to determine suffix)
* @param[in] nvar Number of subvariables (used to determine suffix)
*/
void writeMultiVar( DBfile *fid, const std::string &varname,
const std::vector<std::string> &subVarNames, const std::vector<int> &subVarTypes );

514
analysis/FlowAdaptor.cpp Normal file
View File

@@ -0,0 +1,514 @@
/* Flow adaptor class for multiphase flow methods */
#include "analysis/FlowAdaptor.h"
#include "analysis/distance.h"
#include "analysis/morphology.h"
FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M){
Nx = M.Dm->Nx;
Ny = M.Dm->Ny;
Nz = M.Dm->Nz;
timestep=-1;
timestep_previous=-1;
phi.resize(Nx,Ny,Nz); phi.fill(0); // phase indicator field
phi_t.resize(Nx,Ny,Nz); phi_t.fill(0); // time derivative for the phase indicator field
}
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.01;
double FRACTIONAL_FLOW_EPSILON = 5e-6;
if (M.db->keyExists( "FlowAdaptor" )){
auto flow_db = M.db->getDatabase( "FlowAdaptor" );
MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.01);
FRACTIONAL_FLOW_EPSILON = flow_db->getWithDefault<double>( "fractional_flow_epsilon", 5e-6);
}
int Np = M.Np;
double dA, dB, phi;
double vx,vy,vz;
double mass_a, mass_b, mass_a_global, mass_b_global;
double *Aq_tmp, *Bq_tmp;
double *Vel_x, *Vel_y, *Vel_z, *Phase;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
Phase = new double [Np];
Vel_x = new double [Np];
Vel_y = new double [Np];
Vel_z = new double [Np];
ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double));
ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz;
mass_a = mass_b = 0.0;
double maxSpeed = 0.0;
double localMaxSpeed = 0.0;
/* compute mass change based on weights */
double sum_weights_A = 0.0;
double sum_weights_B = 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++){
int n=M.Map(i,j,k);
//double distance = M.Averages->SDs(i,j,k);
if (!(n<0) ){
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;
vx = Vel_x[n];
vy = Vel_y[n];
vz = Vel_z[n];
double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
double local_weight = (FRACTIONAL_FLOW_EPSILON + local_momentum);
if (phi > 0.0){
sum_weights_A += local_weight*dA;
}
else {
sum_weights_B += local_weight*dB;
}
if ( local_momentum > localMaxSpeed){
localMaxSpeed = local_momentum;
}
}
}
}
}
maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
mass_a_global = M.Dm->Comm.sumReduce(mass_a);
mass_b_global = M.Dm->Comm.sumReduce(mass_b);
double sum_weights_A_global = M.Dm->Comm.sumReduce(sum_weights_A);
double sum_weights_B_global = M.Dm->Comm.sumReduce(sum_weights_B);
sum_weights_A_global /= (FRACTIONAL_FLOW_EPSILON + maxSpeed);
sum_weights_B_global /= (FRACTIONAL_FLOW_EPSILON + maxSpeed);
//double total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
//double total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
/* compute the total mass change */
double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
TOTAL_MASS_CHANGE = 0.1*mass_a_global;
if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
TOTAL_MASS_CHANGE = 0.1*mass_b_global;
double MASS_FACTOR_A = TOTAL_MASS_CHANGE / sum_weights_A_global;
double MASS_FACTOR_B = TOTAL_MASS_CHANGE / sum_weights_B_global;
double LOCAL_MASS_CHANGE = 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++){
int n=M.Map(i,j,k);
if (!(n<0)){
phi = Phase[n];
vx = Vel_x[n];
vy = Vel_y[n];
vz = Vel_z[n];
double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
double local_weight = (FRACTIONAL_FLOW_EPSILON + local_momentum)/(FRACTIONAL_FLOW_EPSILON + maxSpeed);
/* impose ceiling for spurious currents */
//if (local_momentum > maxSpeed) local_momentum = maxSpeed;
if (phi > 0.0){
LOCAL_MASS_CHANGE = MASS_FACTOR_A*local_weight;
Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
//DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
}
else{
LOCAL_MASS_CHANGE = MASS_FACTOR_B*local_weight;
Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
//DebugMassB[n] = LOCAL_MASS_CHANGE;
}
}
}
}
}
if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
return(TOTAL_MASS_CHANGE);
}
void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){
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);
}
double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){
double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.1 );
double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault<double>( "move_interface_factor", 10.0 );
ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) );
/* compute the local derivative of phase indicator field */
double beta = M.beta;
double factor = 0.5/beta;
double total_interface_displacement = 0.0;
double total_interface_sites = 0.0;
for (int n=0; n<Nx*Ny*Nz; n++){
/* compute the distance to the interface */
double value1 = M.Averages->Phi(n);
double dist1 = factor*log((1.0+value1)/(1.0-value1));
double value2 = phi(n);
double dist2 = factor*log((1.0+value2)/(1.0-value2));
phi_t(n) = value2;
if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){
/* time derivative of distance */
double dxdt = 0.125*(dist2-dist1);
/* extrapolate to move the distance further */
double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt;
/* compute the new phase interface */
phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f);
total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt);
total_interface_sites += 1.0;
}
}
ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) );
return total_interface_sites;
}
double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){
const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz);
auto rank = M.rank;
auto Nx = M.Nx; auto Ny = M.Ny; auto Nz = M.Nz;
auto N = Nx*Ny*Nz;
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(M.Dm->Comm,M.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(), M.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 && M.Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
double volume_initial = M.Dm->Comm.sumReduce( count);
double PoreVolume = M.Dm->Volume*M.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;
}
// 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,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm);
M.Dm->Comm.barrier();
// only operate on component "0"ScaLBL_ColorModel &M,
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 = M.Dm->Comm.sumReduce( count);
second_biggest = M.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 (M.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;
}
}
}
}
}
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*M.Dm);
double temp,value;
double factor=0.5/M.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 && M.Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
// erase the original object
phase(i,j,k) = -1.0;
}
}
}
}
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(M.Averages->SDs,phase_distance,phase_id,M.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,*M.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 (M.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*M.beta*d))/(1.f+exp(-2.f*M.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 && M.Averages->SDs(i,j,k) > 0.f){
count+=1.f;
}
}
}
}
double volume_final= M.Dm->Comm.sumReduce( count);
delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial);
if (rank == 0) printf(" new saturation = %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs)));
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double));
// 7. Re-initialize phase field and density
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);
auto BoundaryCondition = M.BoundaryCondition;
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (M.Dm->kproc()==0){
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2);
}
if (M.Dm->kproc() == M.nprocz-1){
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
return delta_volume;
}
double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){
srand(time(NULL));
auto rank = M.rank;
auto Np = M.Np;
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, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
for (int n=0; n < M.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=M.ScaLBL_Comm->FirstInterior(); n < M.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= M.Dm->Comm.sumReduce( count);
mass_loss= M.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(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
return(mass_loss);
}

32
analysis/FlowAdaptor.h Normal file
View File

@@ -0,0 +1,32 @@
/* Flow adaptor class for multiphase flow methods */
#ifndef ScaLBL_FlowAdaptor_INC
#define ScaLBL_FlowAdaptor_INC
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "models/ColorModel.h"
class FlowAdaptor{
public:
FlowAdaptor(ScaLBL_ColorModel &M);
~FlowAdaptor();
double MoveInterface(ScaLBL_ColorModel &M);
double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;
DoubleArray phi_t;
private:
int Nx, Ny, Nz;
int timestep;
int timestep_previous;
};
#endif

26
analysis/analysis.h
View File

@@ -22,8 +22,8 @@ typedef Array<BlobIDType> BlobIDArray;
* @param[in] SignDist SignDist
* @param[in] vF vF
* @param[in] vS vS
* @param[in] S S
* @param[out] LocalBlobID The ids of the blobs
* @param[in] periodic Optional value
* @return Returns the number of blobs
*/
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
@@ -48,12 +48,13 @@ int ComputeLocalPhaseComponent( const IntArray &PhaseID, int &VALUE, IntArray &C
* @param[in] nx Number of elements in the x-direction
* @param[in] ny Number of elements in the y-direction
* @param[in] nz Number of elements in the z-direction
* @param[in] rank_info MPI communication info
* @param[in] Phase Phase
* @param[in] SignDist SignDist
* @param[in] vF vF
* @param[in] vS vS
* @param[in] S S
* @param[out] LocalBlobID The ids of the blobs
* @param[out] GlobalBlobID The ids of the blobs
* @param[in] comm MPI communicator
* @return Returns the number of blobs
*/
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
@@ -68,10 +69,11 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf
* @param[in] nx Number of elements in the x-direction
* @param[in] ny Number of elements in the y-direction
* @param[in] nz Number of elements in the z-direction
* @param[in] rank_in MPI communication info
* @param[in] rank_info MPI communication info
* @param[in] PhaseID Array that identifies the phases
* @param[in] VALUE Identifier for the phase to decompose
* @param[out] GlobalBlobID The ids of the blobs for the phase
* @param[in] comm The communicator to use
* @return Return the number of components in the specified phase
*/
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info,
@@ -82,10 +84,8 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
* @brief Reorder the blobs
* @details Reorder the blobs based on the number of cells they contain
* largest first.
* @param[in] nx Number of elements in the x-direction
* @param[in] ny Number of elements in the y-direction
* @param[in] nz Number of elements in the z-direction
* @param[in/out] ID The ids of the blobs
* @param[in,out] ID The ids of the blobs
* @param[in] comm MPI communicator
*/
void ReorderBlobIDs( BlobIDArray& ID, const Utilities::MPI& comm );
@@ -117,8 +117,12 @@ struct ID_map_struct {
* @details This functions computes the map of blob ids between iterations
* @return Returns the map of the blob ids. Each final blob may have no source
* ids, one parent, or multiple parents. Each src id may be a parent for multiple blobs.
* @param[in] nx Number of elements in the x-direction
* @param[in] ny Number of elements in the y-direction
* @param[in] nz Number of elements in the z-direction
* @param[in] ID1 The blob ids at the first timestep
* @param[in] ID2 The blob ids at the second timestep
* @param[in] comm The communicator to use
*/
ID_map_struct computeIDMap( int nx, int ny, int nz, const BlobIDArray& ID1, const BlobIDArray& ID2, const Utilities::MPI& comm );
@@ -127,7 +131,7 @@ ID_map_struct computeIDMap( int nx, int ny, int nz, const BlobIDArray& ID1, cons
* @brief Compute the new global ids based on the map
* @details This functions computes the time-consistent global ids for the
* current global id index
* @param[in/out] map The timestep mapping for the ids
* @param[in,out] map The timestep mapping for the ids
* @param[in] id_max The globally largest id used previously
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
*/
@@ -139,9 +143,9 @@ void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>&
* @details This functions computes the map of blob ids between iterations.
* Note: we also update the map to reflect the new ids
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
* @param[in/out] IDs The blob ids to renumber
* @param[in,out] IDs The blob ids to renumber
*/
void renumberIDs( const std::vector<BlobIDType>& new_id_list, BlobIDArray& IDs );
void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDArray& IDs );
/*!

2
analysis/distance.h
View File

@@ -24,6 +24,7 @@ inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.
* @param[in] ID Segmentation id
* @param[in] Dm Domain information
* @param[in] periodic Directions that are periodic
* @param[in] dx Cell size
*/
template<class TYPE>
void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
@@ -36,6 +37,7 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
* @param[in] ID Domain id
* @param[in] Dm Domain information
* @param[in] periodic Directions that are periodic
* @param[in] dx Cell size
*/
void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );

3
analysis/filters.h
View File

@@ -25,7 +25,10 @@ void Med3D( const Array<float> &Input, Array<float> &Output );
* @details This routine performs a non-linear local means filter
* @param[in] Input Input image
* @param[in] Mean Mean value
* @param[in] Distance Distance
* @param[out] Output Output image
* @param[in] d
* @param[in] h
*/
int NLM3D( const Array<float> &Input, Array<float> &Mean,
const Array<float> &Distance, Array<float> &Output, const int d, const float h);

3
analysis/imfilter.h
View File

@@ -89,7 +89,7 @@ Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<Array<TY
* multidimensional filter H. The result B has the same size and class as A.
* This version works with separable filters and is more efficient than a single filter.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
* @param[in] Nh The filter size
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
@@ -114,6 +114,7 @@ Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh
* This version works with separable filters and is more efficient than a single filter.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
* @param[in] Nh The filter size
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X

15
analysis/runAnalysis.cpp
View File

@@ -624,8 +624,8 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
d_N[1] = Dm->Ny;
d_N[2] = Dm->Nz;
d_restart_interval = db->getScalar<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_restart_interval = db->getWithDefault<int>( "restart_interval", 100000 );
d_analysis_interval = db->getWithDefault<int>( "analysis_interval", 1000 );
d_subphase_analysis_interval = INT_MAX;
d_visualization_interval = INT_MAX;
d_blobid_interval = INT_MAX;
@@ -639,7 +639,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
}
auto restart_file = db->getScalar<std::string>( "restart_file" );
auto restart_file = db->getWithDefault<std::string>( "restart_file", "Restart");
d_restartFile = restart_file + "." + rankString;
@@ -765,8 +765,8 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
d_N[1] = ColorModel.Dm->Ny;
d_N[2] = ColorModel.Dm->Nz;
d_restart_interval = db->getScalar<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_restart_interval = db->getWithDefault<int>( "restart_interval", 100000 );
d_analysis_interval = db->getWithDefault<int>( "analysis_interval", 1000 );
d_subphase_analysis_interval = INT_MAX;
d_visualization_interval = INT_MAX;
d_blobid_interval = INT_MAX;
@@ -780,10 +780,9 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
}
auto restart_file = db->getScalar<std::string>( "restart_file" );
auto restart_file = db->getWithDefault<std::string>( "restart_file", "Restart");
d_restartFile = restart_file + "." + rankString;
d_rank = d_comm.getRank();
writeIDMap( ID_map_struct(), 0, id_map_filename );
// Initialize IO for silo

2
common/Array.h
View File

@@ -642,7 +642,7 @@ public: // Math operations
void cat( const Array &x, int dim = 0 );
//! Initialize the array with random values (defined from the function table)
void rand();
//void rand();
//! Return true if NaNs are present
bool NaNs() const;

3
common/Array.hpp
View File

@@ -1292,11 +1292,12 @@ TYPE Array<TYPE, FUN, Allocator>::interp( const double *x ) const
/********************************************************
* Math operations (should call the Math class) *
********************************************************/
template<class TYPE, class FUN, class Allocator>
/*template<class TYPE, class FUN, class Allocator>
void Array<TYPE, FUN, Allocator>::rand()
{
FUN::rand( *this );
}
*/
template<class TYPE, class FUN, class Allocator>
Array<TYPE, FUN, Allocator> &
Array<TYPE, FUN, Allocator>::operator+=( const Array<TYPE, FUN, Allocator> &rhs )

1
common/Communication.h
View File

@@ -51,6 +51,7 @@ class fillHalo
public:
/*!
* @brief Default constructor
* @param[in] comm Communicator to use
* @param[in] info Rank and neighbor rank info
* @param[in] n Number of local cells
* @param[in] ng Number of ghost cells

3
common/Domain.cpp
View File

@@ -629,7 +629,8 @@ void Domain::ComputePorosity(){
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocx()*nprocy()*nprocz());
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
if (BoundaryCondition > 0 && BoundaryCondition !=5)
iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-inlet_layers_z - outlet_layers_z));
//.........................................................
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
for (int j=1;j<Ny-1;j++){

4
common/FunctionTable.cpp
View File

@@ -4,7 +4,7 @@
/********************************************************
* Random number generation *
********************************************************/
template<> char genRand<char>()
/*template<> char genRand<char>()
{
static std::random_device rd;
static std::mt19937 gen( rd() );
@@ -88,7 +88,7 @@ template<> long double genRand<long double>()
static std::uniform_real_distribution<double> dis;
return dis( gen );
}
*/
/********************************************************
* axpy *

7
common/FunctionTable.hpp
View File

@@ -7,21 +7,20 @@
#include <algorithm>
#include <cstring>
#include <limits>
#include <random>
//#include <random>
/********************************************************
* Random number initialization *
********************************************************/
template<class TYPE> TYPE genRand();
/*template<class TYPE> TYPE genRand();
template<class TYPE, class FUN>
inline void FunctionTable::rand( Array<TYPE, FUN> &x )
{
for ( size_t i = 0; i < x.length(); i++ )
x( i ) = genRand<TYPE>();
}
*/
/********************************************************
* Reduction *

2
common/MPI.h
View File

@@ -821,7 +821,7 @@ public: // Member functions
* @return Output array for allGather
*/
template<class type>
std::vector<type> allGather( const std::vector<type> &x_in ) const;
std::vector<type> allGather( const std::vector<type> &x ) const;
/*!

2
common/ScaLBL.h
View File

@@ -132,6 +132,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish);
//maybe deprecated
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np);

1
common/Utilities.h
View File

@@ -30,6 +30,7 @@ using StackTrace::Utilities::sleep_s;
* \details This routine will peform the default startup sequence
* \param argc argc from main
* \param argv argv from main
* \param multiple Intialize mpi with MPI_THREAD_MULTIPLE support?
*/
void startup( int argc, char **argv, bool multiple=true );

81
cpu/Poisson.cpp
View File

@@ -235,6 +235,87 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
}
}
extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish){
int n,nn,ijk;
double psi;//electric potential
double rho_e;//local charge density
// neighbors of electric potential psi
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double psi_Laplacian;
double residual_error;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
ijk=Map[n];
psi = Psi[ijk];
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Psi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Psi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Psi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Psi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Psi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Psi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Psi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Psi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Psi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Psi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Psi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Psi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Psi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Psi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Psi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Psi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Psi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Psi[nn]; // get neighbor for phi - 18
psi_Laplacian = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*psi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*psi));//Laplacian of electric potential
residual_error = psi_Laplacian+rho_e/epsilon_LB;
ResidualError[n] = residual_error;
}
}
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np){
//

2
docs/source/conf.py
View File

@@ -18,7 +18,7 @@ import sys
# -- Project information -----------------------------------------------------
project = 'LBPM'
copyright = '2021, James E McClure'
copyright = '2021'
author = 'James E McClure'
# The full version, including alpha/beta/rc tags

7
docs/source/examples/running.rst
View File

@@ -1,9 +1,10 @@
============
Running LBPM
============
==============
LBPM examples
==============
There are two main components to running LBPM simulators.
First is understanding how to launch MPI tasks on your system,
which depends on the particular implementation of MPI that you are using,
as well as other details of the local configuration. The second component is
understanding the LBPM input file structure.

167
docs/source/userGuide/models/color/index.rst
View File

@@ -66,7 +66,172 @@ Model Formulation
****************************
Two LBEs are constructed to model the mass transport, incorporating the anti-diffusion
.. math::
:nowrap:
$$
A_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) = w_q N_a \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
+ \beta \frac{N_b}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big] \;
$$
.. math::
:nowrap:
$$
B_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) =
w_q N_b \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
- \beta \frac{N_a}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big]\;,
$$
The number density for each fluid is obtained from the sum of the mass transport distributions
.. math::
:nowrap:
$$
N_a = \sum_q A_q\;, \quad N_b = \sum_q B_q\;
$$
The phase indicator field is then defined as
.. math::
:nowrap:
$$
\phi = \frac{N_a-N_b}{N_a+N_b}
$$
The fluid density and kinematic viscosity are determined based on linear interpolation
.. math::
:nowrap:
$$
\rho_0 = \frac{(1+\phi) \rho_n}{2}+ \frac{(1-\phi) \rho_w}{2} \;,
$$
.. math::
:nowrap:
$$
s_\nu = \frac{(1+\phi)}{2\tau_n} +\frac{(1-\phi)}{2\tau_w} \;,
$$
where
.. math::
:nowrap:
$$
\nu_w = \frac{1}{3}\Big(\tau_w - \frac{1}{2} \Big) \;, \quad
\nu_n = \frac{1}{3}\Big(\tau_n - \frac{1}{2} \Big) \;.
$$
These values are then used to model the momentum transport.
The LBE governing momentum transport is defined based on a MRT relaxation process with additional
terms to account for the interfacial stresses
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = \sum^{Q-1}_{k=0} M^{-1}_{qk} \lambda_{k} (m_k^{eq}-m_k) + w_q \bm{\xi}_q \cdot \frac{\bm{F}}{c_s^2} \;,
$$
Where :math:`\bm{F}` is an external body force and :math:`c_s^2 = 1/3` is the speed of sound for the LB model.
The moments are linearly indepdendent:
.. math::
:nowrap:
$$
m_k = \sum_{q=0}^{18} M_{qk} f_q\;.
$$
The relaxation parameters are determined from the relaxation time:
.. math::
:nowrap:
$$
\lambda_1 = \lambda_2= \lambda_9 = \lambda_{10}= \lambda_{11}= \lambda_{12}= \lambda_{13}= \lambda_{14}= \lambda_{15} = s_\nu \;,
$$
.. math::
:nowrap:
$$
\lambda_{4}= \lambda_{6}= \lambda_{8} = \lambda_{16} = \lambda_{17} = \lambda_{18}= \frac{8(2-s_\nu)}{8-s_\nu} \;,
$$
The non-zero equilibrium moments are defined as
.. math::
:nowrap:
$$
m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\
$$
.. math::
:nowrap:
$$
m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;,
$$
where the color gradient is determined from the phase indicator field
.. math::
:nowrap:
$$
\textbf{C}=\nabla \phi\;.
$$
and the unit normal vector is
.. math::
:nowrap:
$$
\bm{n} = \frac{\textbf{C}}{|\textbf{C}|}\;.
$$
****************************
Boundary Conditions
****************************

211
docs/source/userGuide/models/mrt/mrt.rst
View File

@@ -1,6 +1,209 @@
=============================================
###############################################################################
MRT model
=============================================
###############################################################################
The LBPM single fluid model is implemented by combining a multi-relaxation time (MRT) D3Q19
lattice Boltzmann equation (LBE) to solve for the momentum transport, recovering the Navier-Stokes
equations to second order based on the Chapman-Enskog expansion. The MRT model is used to assess the
permeability of digital rock images in either the Darcy or non-Darcy flow regimes.
A typical command to launch the LBPM color simulator is as follows
```
mpirun -np $NUMPROCS lbpm_permeability_simulator input.db
```
where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is
the name of the input database that provides the simulation parameters.
Note that the specific syntax to launch MPI tasks may vary depending on your system.
For additional details please refer to your local system documentation.
***************************
Model parameters
***************************
The essential model parameters for the color model are
- ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5`
The kinematic viscosity is given by
.. math::
:nowrap:
$$
\nu = \frac{1}{3} \Big( \tau - \frac 12 \Big)
$$
****************************
Model Formulation
****************************
The LBE governing momentum transport is defined based on a MRT relaxation based on the D3Q19 discrete
velocity set, which determines the values :math:`\bm{\xi}_q`
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) = \sum^{Q-1}_{k=0} M^{-1}_{qk} \lambda_{k} (m_k^{eq}-m_k) + w_q \bm{\xi}_q \cdot \frac{\bm{F}}{c_s^2} \;,
$$
Where :math:`\bm{F}` is an external body force and :math:`c_s^2 = 1/3` is the speed of sound for the LB model.
The moments are linearly indepdendent functions of the distributions:
.. math::
:nowrap:
$$
m_k = \sum_{q=0}^{18} M_{qk} f_q\;.
$$
The non-zero equilibrium moments are
.. math::
:nowrap:
$$
m_1^{eq} = 19\frac{j_x^2+j_y^2+j_z^2}{\rho} - 11\rho \;,
$$
.. math::
:nowrap:
$$
m_2^{eq} = 3\rho - \frac{11}{2} \frac{j_x^2+j_y^2+j_z^2}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_4^{eq} = -\frac 2 3 j_x \;,
$$
.. math::
:nowrap:
$$
m_6^{eq} = -\frac 2 3 j_y \;,
$$
.. math::
:nowrap:
$$
m_8^{eq} = -\frac 2 3 j_z \;,
$$
.. math::
:nowrap:
$$
m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho}\;,
$$
.. math::
:nowrap:
$$
m_{10}^{eq} = -\frac{2j_x^2-j_y^2-j_z^2)}{2\rho} \;,
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{12}^{eq} = -\frac{j_y^2-j_z^2}{2\rho} \;,
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = \frac{j_x j_y}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = \frac{j_y j_z}{\rho} \;,
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = \frac{j_x j_z}{\rho} \;,
$$
The relaxation parameters are determined based on the relaxation time :math:`\tau`
.. math::
:nowrap:
$$
\lambda_1 = \lambda_2= \lambda_9 = \lambda_{10}= \lambda_{11}= \lambda_{12}= \lambda_{13}= \lambda_{14}= \lambda_{15} = s_\nu = \frac{1}{\tau} \;,
$$
.. math::
:nowrap:
$$
\lambda_{4}= \lambda_{6}= \lambda_{8} = \lambda_{16} = \lambda_{17} = \lambda_{18}= \frac{8(2-s_\nu)}{8-s_\nu} \;,
$$
****************************
Example Input File
****************************
****************************
Boundary Conditions
****************************
The following external boundary conditions are supported by ``lbpm_permeability_simulator``
and can be set by setting the ``BC`` key values in the ``Domain`` section of the
input file database
- ``BC = 0`` -- fully periodic boundary conditions
- ``BC = 3`` -- constant pressure boundary condition
- ``BC = 4`` -- constant volumetric flux boundary condition
For ``BC = 0`` any mass that exits on one side of the domain will re-enter at the other
side. If the pore-structure for the image is tight, the mismatch between the inlet and
outlet can artificially reduce the permeability of the sample due to the blockage of
flow pathways at the boundary. LBPM includes an internal utility that will reduce the impact
of the boundary mismatch by eroding the solid labels within the inlet and outlet layers
(https://doi.org/10.1007/s10596-020-10028-9) to create a mixing layer.
The number mixing layers to use can be set using the key values in the ``Domain`` section
of the input database
- ``InletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the inlet
- ``OUtletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the outlet
For the other boundary conditions a thin reservoir of fluid (default ``3`` voxels)
is established at either side of the domain. The inlet is defined as the boundary face
where ``z = 0`` and the outlet is the boundary face where ``z = nprocz*nz``. By default a
reservoir of fluid A is established at the inlet and a reservoir of fluid B is established at
the outlet, each with a default thickness of three voxels. To over-ride the default label at
the inlet or outlet, the ``Domain`` section of the database may specify the following key values
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet
The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the
solution of the Navier-Stokes equations.

33
doxygen/Doxyfile.in
View File

@@ -379,12 +379,6 @@ MAX_INITIALIZER_LINES = 30
SHOW_USED_FILES = YES
# If the sources in your project are distributed over multiple directories
# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy
# in the documentation. The default is NO.
SHOW_DIRECTORIES = YES
# Set the SHOW_FILES tag to NO to disable the generation of the Files page.
# This will remove the Files entry from the Quick Index and from the
# Folder Tree View (if specified). The default is YES.
@@ -467,7 +461,7 @@ INPUT_ENCODING = UTF-8
# *.hxx *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.dox *.py
# *.f90 *.f *.for *.vhd *.vhdl
FILE_PATTERNS = *.h *.hh *.cu
FILE_PATTERNS = *.h *.hh
# The RECURSIVE tag can be used to turn specify whether or not subdirectories
# should be searched for input files as well. Possible values are YES and NO.
@@ -728,12 +722,6 @@ HTML_COLORSTYLE_GAMMA = 80
HTML_TIMESTAMP = YES
# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
# files or namespaces will be aligned in HTML using tables. If set to
# NO a bullet list will be used.
HTML_ALIGN_MEMBERS = YES
# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML
# documentation will contain sections that can be hidden and shown after the
# page has loaded. For this to work a browser that supports
@@ -850,7 +838,7 @@ TREEVIEW_WIDTH = 250
# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
# generate Latex output.
GENERATE_LATEX = YES
GENERATE_LATEX = NO
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
# If a relative path is entered the value of OUTPUT_DIRECTORY will be
@@ -1005,18 +993,6 @@ GENERATE_XML = NO
XML_OUTPUT = xml
# The XML_SCHEMA tag can be used to specify an XML schema,
# which can be used by a validating XML parser to check the
# syntax of the XML files.
XML_SCHEMA =
# The XML_DTD tag can be used to specify an XML DTD,
# which can be used by a validating XML parser to check the
# syntax of the XML files.
XML_DTD =
# If the XML_PROGRAMLISTING tag is set to YES Doxygen will
# dump the program listings (including syntax highlighting
# and cross-referencing information) to the XML output. Note that
@@ -1174,11 +1150,6 @@ ALLEXTERNALS = NO
EXTERNAL_GROUPS = YES
# The PERL_PATH should be the absolute path and name of the perl script
# interpreter (i.e. the result of `which perl').
PERL_PATH = /usr/bin/perl
#---------------------------------------------------------------------------
# Configuration options related to the dot tool
#---------------------------------------------------------------------------

2
doxygen/DoxygenMainpage.h
View File

@@ -4,8 +4,6 @@
*
* - \ref IO "IO routines"
* - \ref Utilities "Utility routines"
* - \ref silo "Access to silo routines"
* - \ref netcdf "Access to netcdf routines"
*
* \author James McClure
*/

6
example/Bubble/input.db
View File

@@ -47,4 +47,8 @@ Analysis {
}
Visualization {
}
}
FlowAdaptor {
}

2
example/MicroModel/BlobSize.in
View File

@@ -1,2 +0,0 @@
1
32 32 32

6
example/MicroModel/Color.in
View File

@@ -1,6 +0,0 @@
0.7
1.0e-2 0.95 -1.0
0.7
0.0 0.0 0.0
0 1 1.04 0.96
900000 20000 1e-5

144
example/MicroModel/DiscPack.in
View File

@@ -1,144 +0,0 @@
33.13248649 29.93643243 15.41628378
75.07072973 26.66563068 16.99835135
118.9166622 26.51018919 13.54747297
155.9032432 26.54839365 18.37214865
194.4647297 25.64848784 13.23047838
236.4466216 27.04415726 18.86951351
276.9660811 27.18572081 18.66178378
316.6336486 27.16016042 17.04760811
360.6031081 29.63283514 14.62525676
400.4937838 26.33180541 15.91213514
436.6832432 29.38145946 15.30281081
479.5852703 27.09887028 18.76375676
15.17307297 69.60128378 14.82944595
53.0137973 64.87014865 16.21706757
94.52747297 68.84901351 17.6057027
135.5813919 65.70918919 16.92474324
173.2494595 69.96463514 14.05135135
215.6589189 68.16136486 17.94732432
259.9767568 70.38693243 14.80982432
297.0702703 67.83155405 14.68744595
336.5431081 67.57271622 17.27831081
378.3740541 67.72916216 18.13277027
419.2805405 66.79117568 17.42808108
455.1844595 72.06137838 13.27178784
33.72422973 108.8409189 17.24389189
74.43232432 108.0137838 18.75
114.1253784 106.4523919 17.1622027
153.882027 112.3773514 13.02102703
199.0522973 111.1353243 15.15556757
241.1697297 107.2128378 13.80640541
277.315 108.4926486 18.42555405
320.4716216 106.1688919 15.76947297
354.8572973 105.7547973 13.44669865
400.797973 110.3642297 16.28766216
439.5955405 107.9637838 18.33110811
482.8833784 107.3503784 13.76293243
10.30476351 149.3496351 15.64486486
55.28555405 147.8185946 16.92521622
95.24567568 149.0092432 17.58089189
131.8110405 150.626527 13.8865
179.3263514 146.4626757 14.22206757
216.0708108 152.3747568 14.82614865
261.7548649 147.6711757 13.01481622
299.3959459 147.4904324 13.76047297
339.8782432 149.7220541 16.66758108
373.3464865 152.5632297 13.45467297
417.9712162 147.855973 17.36840541
458.4025676 144.5771622 14.61705405
33.32618919 191.2175676 12.4216527
74.30910811 189.2091892 18.52312162
114.6439054 190.2966216 16.95763514
158.2144595 191.8867568 15.55116216
194.6536486 189.0091892 15.71814865
237.9637838 189.242027 17.27986486
271.4318919 186.4898649 13.0739527
317.7631081 188.4841892 16.13393243
358.1116216 191.2745946 14.83137838
398.7190541 189.4286486 18.13921622
439.5593243 189.8759459 16.8382973
482.1985135 194.4328378 12.19701081
13.90195446 229.4845946 18.36889189
55.69509459 226.9613514 16.02147297
96.98522973 232.4032432 14.46966216
137.4593243 227.1902703 13.4368527
174.2690541 228.4505405 16.73159459
216.2844595 229.237027 18.25581081
259.6358108 233.6717568 12.89051757
298.1172973 229.1837838 15.79105405
338.9263514 229.9989189 16.83474324
378.4617568 230.0036486 18.49422973
418.9406757 228.4566216 13.19394865
459.3066216 230.6732432 17.26921622
35.56972973 273.1122973 15.46828378
72.91724324 272.0410811 16.95751351
114.3243243 273.7995946 15.17266216
158.1804054 267.382973 13.8992027
194.7917568 269.7763514 17.70321622
236.57 271.7548649 12.18126081
272.412973 266.6922973 13.18252973
317.5945946 270.2414865 18.8062027
359.0810811 268.9272973 17.10451351
401.6145946 268.1164865 13.57114865
434.9897297 265.9721622 12.38897297
479.9452703 270.4847297 12.98088378
13.27824108 310.6754054 18.66364865
53.25386486 312.0385135 14.46228378
92.49813514 308.3451351 15.20474324
136.1003514 309.7871622 17.76025676
174.5641892 311.4295946 13.06244865
219.4805405 310.8893243 14.40940541
258.0839189 310.0159459 17.41775676
297.6322973 310.9914865 18.51806757
337.9448649 310.7555405 17.73358108
374.9710811 312.2205405 14.25995946
422.8451351 314.0112162 14.40621622
459.8067568 311.2154054 13.47336486
34.99885135 351.7381081 17.47081081
75.33945946 351.677027 15.52572973
114.8586486 351.3459459 18.89374324
155.7139189 350.8706757 18.07736486
191.6158108 352.282973 12.90032838
234.2532432 350.5644595 14.65478378
275.127027 350.5708108 16.79345946
317.652973 352.3060811 16.43624324
363.8067568 347.8472973 12.7394473
401.5713514 349.2537838 14.74306757
436.0778378 349.6571622 14.30860811
482.565 353.0252703 15.15274324
10.92297703 391.8589189 15.69475676
52.44290541 394.7612162 15.2132027
93.64185135 390.0597297 15.0062027
139.127527 389.3374324 14.73910811
175.5394595 391.7282432 18.06841892
216.4606757 392.4894595 17.91210811
251.7971622 390.4447297 12.52764324
298.7958108 393.1183784 16.58391892
336.7648649 392.7705405 16.06710811
379.9774324 392.8135135 16.87058108
417.1693243 392.1031081 16.41077027
459.4625676 391.9014865 18.88939189
32.60527027 431.387973 13.79572973
73.08481081 428.7002703 14.2487973
114.5770541 432.5933784 17.63783784
156.0283784 432.7117568 17.82533784
198.9813514 432.5832432 14.1537027
237.5352703 434.3167568 15.54622973
277.7677027 433.0071622 15.02437838
316.0105405 434.197973 16.34566216
359.23 437.6654054 13.07579189
401.887027 430.8952703 15.38594595
442.0587838 434.6848649 12.87081892
481.2160811 438.0641892 13.1349027
13.43983284 473.0802703 17.79305405
58.56352703 468.0102703 12.698
95.47004054 472.9767568 16.32575676
134.9401622 472.807027 18.70305405
173.3897297 469.9112162 15.53236486
216.2163514 473.1051351 18.61358108
256.3348649 474.0559459 14.83377027
294.8790541 471.1768919 14.53377027
339.6306757 472.1563514 16.27194595
377.4156757 471.6939189 17.32056757
423.6267568 474.1678378 12.76715541
459.5241892 471.4267568 16.65775676

3
example/MicroModel/Domain.in
View File

@@ -1,3 +0,0 @@
1 1 1
10 400 400
20.0 500.0 500.0

5466
example/MicroModel/FullMicromodel.discs Normal file
View File

File diff suppressed because it is too large Load Diff

47
example/MicroModel/GenerateCases.sh Executable file
View File

@@ -0,0 +1,47 @@
#!/bin/bash
tau1=1.18
tau2=0.7
alpha=0.005
Q="1.179757e-08 1.179757e-07 1.179757e-06"
#Cases for drainage
DrainWet="0.79 0.47 0.0"
# Cases for imbibition
ImbWet="0.92 0.47"
for q in $Q; do
echo $q;
flux=$(echo $q | sed 's/1.179757e-08/0.002/g')
flux=$(echo $flux | sed 's/1.179757e-07/0.02/g')
flux=$(echo $flux | sed 's/1.179757e-06/0.2/g')
for i in $DrainWet; do
NAME="Juanes_drain_Q"$flux"_wet"$i
echo $NAME
mkdir $NAME
echo "$tau1 $tau2" > $NAME/Color.in
echo "$alpha 0.95 $i" >> $NAME/Color.in
echo "0.0" >> $NAME/Color.in
echo "0.0 0.0 0.0" >> $NAME/Color.in
echo "0 4 $q 0.0" >> $NAME/Color.in
echo "5000000 25000 1e-5" >> $NAME/Color.in
echo "1000" >> $NAME/Color.in
done
for i in $ImbWet; do
NAME="Juanes_imb_Q"$flux"_wet"$i
echo $NAME
mkdir $NAME
echo "$tau1 $tau2" > $NAME/Color.in
echo "$alpha 0.95 $i" >> $NAME/Color.in
echo "0.0" >> $NAME/Color.in
echo "0.0 0.0 0.0" >> $NAME/Color.in
echo "0 4 $q 0.0" >> $NAME/Color.in
echo "5000000 25000 1e-5" >> $NAME/Color.in
echo "1000" >> $NAME/Color.in
done
done

21
example/MicroModel/ReadDiscPack.R Normal file
View File

@@ -0,0 +1,21 @@
require("ggplot2")
Discs<-read.csv("FullMicromodel.discs",head=FALSE,sep=" ")
colnames(Discs)<-c("cx","cy","radius")
L=0.45
# Extract some subset from the interior of the discs
SubDiscs<-subset(Discs,Discs$cy>0.9-L)
SubDiscs<-subset(SubDiscs,SubDiscs$cy<0.9+L)
SubDiscs<-subset(SubDiscs,SubDiscs$cx>0.9-L)
SubDiscs<-subset(SubDiscs,SubDiscs$cx<0.9+L)
SubDiscs$cx<-SubDiscs$cx-0.9+L
SubDiscs$cy<-SubDiscs$cy-0.9+L
write.table(SubDiscs,file="DiscPack.in",quote=FALSE,row.names=FALSE,col.names=FALSE,sep=" ")
ShowPlot<-ggplot(SubDiscs)+geom_circle(aes(x0=cx,y0=cy,r=radius))

6
example/Piston/input.db
View File

@@ -37,4 +37,8 @@ Analysis {
Visualization {
}
}
FlowAdaptor {
}

6
example/Plates/input.db
View File

@@ -36,4 +36,8 @@ Analysis {
}
Visualization {
}
}
FlowAdaptor {
}

617
models/ColorModel.cpp
View File

@@ -198,22 +198,12 @@ void ScaLBL_ColorModel::ReadParams(string filename){
domain_db->putScalar<int>( "BC", BoundaryCondition );
}
else if (protocol == "core flooding"){
if (rank == 0) printf("Using core flooding protocol \n");
if (rank == 0) printf("Using core flooding protocol \n");
if (BoundaryCondition != 4){
BoundaryCondition = 4;
if (rank==0) printf("WARNING: protocol (core flooding) supports only volumetric flux boundary condition \n");
}
domain_db->putScalar<int>( "BC", BoundaryCondition );
if (color_db->keyExists( "capillary_number" )){
double capillary_number = color_db->getScalar<double>( "capillary_number" );
if (rank==0) printf(" set flux to achieve Ca=%f \n", capillary_number);
double MuB = rhoB*(tauB - 0.5)/3.0;
double IFT = 6.0*alpha;
//double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
flux = Dm->Porosity()*nprocx*(Nx-2)*nprocy*(Ny-2)*IFT*capillary_number/MuB;
if (rank==0) printf(" flux=%f \n",flux);
}
color_db->putScalar<double>( "flux", flux );
}
}
@@ -493,6 +483,18 @@ void ScaLBL_ColorModel::Create(){
void ScaLBL_ColorModel::Initialize(){
/* if both capillary number and flux BC are specified */
if (color_db->keyExists( "capillary_number" ) && BoundaryCondition == 4){
double capillary_number = color_db->getScalar<double>( "capillary_number" );
if (rank==0) printf(" set flux to achieve Ca=%f \n", capillary_number);
double MuB = rhoB*(tauB - 0.5)/3.0;
double IFT = 6.0*alpha;
double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
flux = Mask->Porosity()*CrossSectionalArea*(Ny-2)*IFT*capillary_number/MuB;
if (rank==0) printf(" flux=%f \n",flux);
}
color_db->putScalar<double>( "flux", flux );
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
/*
@@ -588,19 +590,24 @@ double ScaLBL_ColorModel::Run(int returntime){
bool Regular = false;
bool RESCALE_FORCE = false;
bool SET_CAPILLARY_NUMBER = false;
bool TRIGGER_FORCE_RESCALE = false;
double tolerance = 0.01;
auto current_db = db->cloneDatabase();
auto flow_db = db->getDatabase( "FlowAdaptor" );
int MIN_STEADY_TIMESTEPS = flow_db->getWithDefault<int>( "min_steady_timesteps", 1000000 );
int MAX_STEADY_TIMESTEPS = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
int RESCALE_FORCE_AFTER_TIMESTEP = MAX_STEADY_TIMESTEPS*2;
int INITIAL_TIMESTEP = timestep;
double capillary_number = 1.0e-5;
double Ca_previous = 0.0;
double minCa = 8.0e-6;
double maxCa = 1.0;
if (color_db->keyExists( "capillary_number" )){
capillary_number = color_db->getScalar<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
maxCa = 2.0*capillary_number;
minCa = 0.8*capillary_number;
}
if (color_db->keyExists( "rescale_force_after_timestep" )){
RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar<int>( "rescale_force_after_timestep" );
@@ -737,8 +744,26 @@ double ScaLBL_ColorModel::Run(int returntime){
isSteady = true;
if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS)
isSteady = true;
if (isSteady && (Ca > maxCa || Ca < minCa) && SET_CAPILLARY_NUMBER ){
/* re-run the point if the actual Ca is too far from the target Ca */
isSteady = false;
RESCALE_FORCE = true;
t1 = std::chrono::system_clock::now();
CURRENT_TIMESTEP = 0;
timestep = INITIAL_TIMESTEP;
TRIGGER_FORCE_RESCALE = true;
if (rank == 0) printf(" Capillary number missed target value = %f (measured value was Ca = %f) \n ",capillary_number, Ca);
}
if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP){
TRIGGER_FORCE_RESCALE = true;
}
if (TRIGGER_FORCE_RESCALE){
RESCALE_FORCE = false;
TRIGGER_FORCE_RESCALE = 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;
@@ -2009,571 +2034,3 @@ void ScaLBL_ColorModel::WriteDebug(){
*/
}
FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M){
Nx = M.Dm->Nx;
Ny = M.Dm->Ny;
Nz = M.Dm->Nz;
timestep=-1;
timestep_previous=-1;
phi.resize(Nx,Ny,Nz); phi.fill(0); // phase indicator field
phi_t.resize(Nx,Ny,Nz); phi_t.fill(0); // time derivative for the phase indicator field
}
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;
if (M.db->keyExists( "FlowAdaptor" )){
auto flow_db = M.db->getDatabase( "FlowAdaptor" );
MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.0005);
}
int Np = M.Np;
double dA, dB, phi;
double vx,vy,vz;
double vax,vay,vaz;
double vbx,vby,vbz;
double vax_global,vay_global,vaz_global;
double vbx_global,vby_global,vbz_global;
double mass_a, mass_b, mass_a_global, mass_b_global;
double *Aq_tmp, *Bq_tmp;
double *Vel_x, *Vel_y, *Vel_z, *Phase;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
Phase = new double [Np];
Vel_x = new double [Np];
Vel_y = new double [Np];
Vel_z = new double [Np];
ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double));
ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz;
/* compute the total momentum */
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++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int n=M.Map(i,j,k);
double distance = M.Averages->SDs(i,j,k);
if (!(n<0) ){
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];
}
double speed = sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz);
if (distance > 3.0 && speed > localMaxSpeed){
localMaxSpeed = speed;
}
}
}
}
}
maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
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);
vay_global = M.Dm->Comm.sumReduce(vay);
vaz_global = M.Dm->Comm.sumReduce(vaz);
vbx_global = M.Dm->Comm.sumReduce(vbx);
vby_global = M.Dm->Comm.sumReduce(vby);
vbz_global = M.Dm->Comm.sumReduce(vbz);
double total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
double total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
/* compute the total mass change */
double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
TOTAL_MASS_CHANGE = 0.1*mass_a_global;
if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
TOTAL_MASS_CHANGE = 0.1*mass_b_global;
double LOCAL_MASS_CHANGE = 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++){
int n=M.Map(i,j,k);
if (!(n<0)){
phi = Phase[n];
vx = Vel_x[n];
vy = Vel_y[n];
vz = Vel_z[n];
double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
/* impose ceiling for spurious currents */
if (local_momentum > maxSpeed) local_momentum = maxSpeed;
if (phi > 0.0){
LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
//DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
}
else{
LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
//DebugMassB[n] = LOCAL_MASS_CHANGE;
}
}
}
}
}
if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
return(TOTAL_MASS_CHANGE);
}
void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){
int Np = M.Np;
double dA, dB;
double *Aq_tmp, *Bq_tmp;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
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];
if (dA > 1.0){
double mass_change = dA - 1.0;
Aq_tmp[n] -= 0.333333333333333*mass_change;
Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
if (dB > 1.0){
double mass_change = dB - 1.0;
Bq_tmp[n] -= 0.333333333333333*mass_change;
Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
}
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];
if (dA > 1.0){
double mass_change = dA - 1.0;
Aq_tmp[n] -= 0.333333333333333*mass_change;
Aq_tmp[n+Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
if (dB > 1.0){
double mass_change = dB - 1.0;
Bq_tmp[n] -= 0.333333333333333*mass_change;
Bq_tmp[n+Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change;
Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change;
}
}
ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
}
double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){
double INTERFACE_CUTOFF = M.color_db->getWithDefault<double>( "move_interface_cutoff", 0.1 );
double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault<double>( "move_interface_factor", 10.0 );
ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) );
/* compute the local derivative of phase indicator field */
double beta = M.beta;
double factor = 0.5/beta;
double total_interface_displacement = 0.0;
double total_interface_sites = 0.0;
for (int n=0; n<Nx*Ny*Nz; n++){
/* compute the distance to the interface */
double value1 = M.Averages->Phi(n);
double dist1 = factor*log((1.0+value1)/(1.0-value1));
double value2 = phi(n);
double dist2 = factor*log((1.0+value2)/(1.0-value2));
phi_t(n) = value2;
if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){
/* time derivative of distance */
double dxdt = 0.125*(dist2-dist1);
/* extrapolate to move the distance further */
double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt;
/* compute the new phase interface */
phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f);
total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt);
total_interface_sites += 1.0;
}
}
ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) );
return total_interface_sites;
}
double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){
const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz);
auto rank = M.rank;
auto Nx = M.Nx; auto Ny = M.Ny; auto Nz = M.Nz;
auto N = Nx*Ny*Nz;
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(M.Dm->Comm,M.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(), M.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 && M.Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
double volume_initial = M.Dm->Comm.sumReduce( count);
double PoreVolume = M.Dm->Volume*M.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;
}
// 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,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm);
M.Dm->Comm.barrier();
// only operate on component "0"ScaLBL_ColorModel &M,
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 = M.Dm->Comm.sumReduce( count);
second_biggest = M.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 (M.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;
}
}
}
}
}
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*M.Dm);
double temp,value;
double factor=0.5/M.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 && M.Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
// erase the original object
phase(i,j,k) = -1.0;
}
}
}
}
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(M.Averages->SDs,phase_distance,phase_id,M.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,*M.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 (M.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*M.beta*d))/(1.f+exp(-2.f*M.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 && M.Averages->SDs(i,j,k) > 0.f){
count+=1.f;
}
}
}
}
double volume_final= M.Dm->Comm.sumReduce( count);
delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial);
if (rank == 0) printf(" new saturation = %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs)));
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double));
// 7. Re-initialize phase field and density
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);
auto BoundaryCondition = M.BoundaryCondition;
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
if (M.Dm->kproc()==0){
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2);
}
if (M.Dm->kproc() == M.nprocz-1){
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
return delta_volume;
}
double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){
srand(time(NULL));
auto rank = M.rank;
auto Np = M.Np;
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, M.Aq, 7*Np*sizeof(double));
ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
for (int n=0; n < M.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=M.ScaLBL_Comm->FirstInterior(); n < M.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= M.Dm->Comm.sumReduce( count);
mass_loss= M.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(M.Aq, Aq_tmp, 7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
return(mass_loss);
}

18
models/ColorModel.h
View File

@@ -10,6 +10,7 @@ Implementation of color lattice boltzmann model
#include <fstream>
#include "common/Communication.h"
#include "analysis/FlowAdaptor.h"
#include "analysis/TwoPhase.h"
#include "analysis/runAnalysis.h"
#include "common/MPI.h"
@@ -93,22 +94,5 @@ private:
double MorphOpenConnected(double target_volume_change);
};
class FlowAdaptor{
public:
FlowAdaptor(ScaLBL_ColorModel &M);
~FlowAdaptor();
double MoveInterface(ScaLBL_ColorModel &M);
double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;
DoubleArray phi_t;
private:
int Nx, Ny, Nz;
int timestep;
int timestep_previous;
};
#endif

28
models/IonModel.cpp
View File

@@ -627,13 +627,13 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
}
}
void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion)
void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion, int ic)
{
double *Ci_host;
Ci_host = new double[N];
double VALUE=0.f;
Mask->ReadFromFile(File_ion[0],File_ion[1],Ci_host);
Mask->ReadFromFile(File_ion[2*ic],File_ion[2*ic+1],Ci_host);
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
@@ -722,20 +722,24 @@ void ScaLBL_IonModel::Initialize(){
*/
if (rank==0) printf ("LB Ion Solver: initializing D3Q7 distributions\n");
if (ion_db->keyExists("IonConcentrationFile")){
//TODO: Need to figure out how to deal with multi-species concentration initialization
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
double *Ci_host;
Ci_host = new double[number_ion_species*Np];
for (size_t ic=0; ic<number_ion_species; ic++){
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
if (File_ion.size()==2*number_ion_species){
double *Ci_host;
Ci_host = new double[number_ion_species*Np];
for (size_t ic=0; ic<number_ion_species; ic++){
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion,ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
comm.barrier();
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
}
delete [] Ci_host;
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
comm.barrier();
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
else{
ERROR("Error: Number of user-input ion concentration files should be equal to number of ion species!\n");
}
delete [] Ci_host;
}
else{
for (size_t ic=0; ic<number_ion_species; ic++){

2
models/IonModel.h
View File

@@ -96,7 +96,7 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *ion_solid);
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion);
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion,int ic);
void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
};
#endif

116
models/PoissonSolver.cpp
View File

@@ -69,6 +69,8 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "tolerance" )){
tolerance = electric_db->getScalar<double>( "tolerance" );
}
//'tolerance_method' can be {"MSE","MSE_max"}
tolerance_method = electric_db->getWithDefault<std::string>( "tolerance_method", "MSE" );
if (electric_db->keyExists( "epsilonR" )){
epsilonR = electric_db->getScalar<double>( "epsilonR" );
}
@@ -122,6 +124,15 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance);
if (rank==0) printf(" LB relaxation tau = %.5g \n", tau);
if (rank==0) printf("***********************************************************************************\n");
if (tolerance_method.compare("MSE")==0){
if (rank==0) printf("LB-Poisson Solver: Use averaged MSE to check solution convergence.\n");
}
else if (tolerance_method.compare("MSE_max")==0){
if (rank==0) printf("LB-Poisson Solver: Use maximum MSE to check solution convergence.\n");
}
else{
if (rank==0) printf("LB-Poisson Solver: tolerance_method=%s cannot be identified!\n",tolerance_method.c_str());
}
switch (BoundaryConditionSolid){
case 1:
@@ -150,6 +161,7 @@ void ScaLBL_Poisson::SetDomain(){
N = Nx*Ny*Nz;
Distance.resize(Nx,Ny,Nz);
Psi_host.resize(Nx,Ny,Nz);
Psi_previous.resize(Nx,Ny,Nz);
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
@@ -338,6 +350,7 @@ void ScaLBL_Poisson::Create(){
ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &ResidualError, sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
@@ -541,12 +554,10 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
timestep=0;
double error = 1.0;
double psi_avg_previous = 0.0;
while (timestep < timestepMax && error > tolerance) {
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
SolvePoissonAAodd(ChargeDensity);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
@@ -559,33 +570,86 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
//************************************************************************/
// Check convergence of steady-state solution
if (timestep==2){
//save electric potential for convergence check
ScaLBL_CopyToHost(Psi_previous.data(),Psi,sizeof(double)*Nx*Ny*Nz);
}
if (timestep%analysis_interval==0){
//ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
double count_loc=0;
double count;
double psi_avg;
double psi_loc=0.f;
if (tolerance_method.compare("MSE")==0){
double count_loc=0;
double count;
double MSE_loc=0.0;
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
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 (Distance(i,j,k) > 0){
MSE_loc += (Psi_host(i,j,k) - Psi_previous(i,j,k))*(Psi_host(i,j,k) - Psi_previous(i,j,k));
count_loc+=1.0;
}
}
}
}
error=Dm->Comm.sumReduce(MSE_loc);
count=Dm->Comm.sumReduce(count_loc);
error /= count;
}
else if (tolerance_method.compare("MSE_max")==0){
vector<double>MSE_loc;
double MSE_loc_max;
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
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 (Distance(i,j,k) > 0){
MSE_loc.push_back((Psi_host(i,j,k) - Psi_previous(i,j,k))*(Psi_host(i,j,k) - Psi_previous(i,j,k)));
}
}
}
}
vector<double>::iterator it_max = max_element(MSE_loc.begin(),MSE_loc.end());
unsigned int idx_max=distance(MSE_loc.begin(),it_max);
MSE_loc_max=MSE_loc[idx_max];
error=Dm->Comm.maxReduce(MSE_loc_max);
}
else{
ERROR("Error: user-specified tolerance_method cannot be identified; check you input database! \n");
}
ScaLBL_CopyToHost(Psi_previous.data(),Psi,sizeof(double)*Nx*Ny*Nz);
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 (Distance(i,j,k) > 0){
psi_loc += Psi_host(i,j,k);
count_loc+=1.0;
}
}
}
}
psi_avg=Dm->Comm.sumReduce( psi_loc);
count=Dm->Comm.sumReduce( count_loc);
psi_avg /= count;
double psi_avg_mag=psi_avg;
if (psi_avg==0.0) psi_avg_mag=1.0;
error = fabs(psi_avg-psi_avg_previous)/fabs(psi_avg_mag);
psi_avg_previous = psi_avg;
//legacy code that tried to use residual to check convergence
//ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior());
//ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,0, ScaLBL_Comm->LastExterior());
//ScaLBL_Comm->Barrier(); comm.barrier();
//vector<double> ResidualError_host(Np);
//double error_loc_max;
////calculate the maximum residual error
//ScaLBL_CopyToHost(&ResidualError_host[0],ResidualError,sizeof(double)*Np);
//vector<double>::iterator it_temp1,it_temp2;
//it_temp1=ResidualError_host.begin();
//advance(it_temp1,ScaLBL_Comm->LastExterior());
//vector<double>::iterator it_max = max_element(ResidualError_host.begin(),it_temp1);
//unsigned int idx_max1 = distance(ResidualError_host.begin(),it_max);
//it_temp1=ResidualError_host.begin();
//it_temp2=ResidualError_host.begin();
//advance(it_temp1,ScaLBL_Comm->FirstInterior());
//advance(it_temp2,ScaLBL_Comm->LastInterior());
//it_max = max_element(it_temp1,it_temp2);
//unsigned int idx_max2 = distance(ResidualError_host.begin(),it_max);
//if (ResidualError_host[idx_max1]>ResidualError_host[idx_max2]){
// error_loc_max=ResidualError_host[idx_max1];
//}
//else{
// error_loc_max=ResidualError_host[idx_max2];
//}
//error = Dm->Comm.maxReduce(error_loc_max);
}
}
if(WriteLog==true){

4
models/PoissonSolver.h
View File

@@ -10,6 +10,7 @@
#include <stdexcept>
#include <fstream>
#include <cmath>
#include <iterator>
#include "common/ScaLBL.h"
#include "common/Communication.h"
@@ -48,6 +49,7 @@ public:
int BoundaryConditionSolid;
double tau;
double tolerance;
std::string tolerance_method;
double k2_inv;
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;
double Vin, Vout;
@@ -78,6 +80,7 @@ public:
IntArray Map;
DoubleArray Distance;
DoubleArray Psi_host;
DoubleArray Psi_previous;
int *NeighborList;
int *dvcMap;
//signed char *dvcID;
@@ -85,6 +88,7 @@ public:
double *Psi;
double *ElectricField;
double *ChargeDensityDummy;// for debugging
double *ResidualError;
private:
Utilities::MPI comm;

1
tests/CMakeLists.txt
View File

@@ -56,6 +56,7 @@ ADD_LBPM_TEST( TestTorus )
ADD_LBPM_TEST( TestTorusEvolve )
ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestFlowAdaptor )
ADD_LBPM_TEST( TestMap )
#ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad )

254
tests/TestFlowAdaptor.cpp Normal file
View File

@@ -0,0 +1,254 @@
/* Test the flow adaptor class with color model */
#include <exception>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include "common/Utilities.h"
#include "models/ColorModel.h"
inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
// initialize a bubble
int i,j,k,n;
int rank = ColorModel.Dm->rank();
int nprocx = ColorModel.Dm->nprocx();
int nprocy = ColorModel.Dm->nprocy();
int nprocz = ColorModel.Dm->nprocz();
int Nx = ColorModel.Dm->Nx;
int Ny = ColorModel.Dm->Ny;
int Nz = ColorModel.Dm->Nz;
if (rank == 0) cout << "Setting up bubble..." << endl;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
ColorModel.Averages->SDs(i,j,k) = 100.f;
}
}
}
// Initialize the bubble
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;
double iglobal= double(i+(Nx-2)*ColorModel.Dm->iproc())-double((Nx-2)*nprocx)*0.5;
double jglobal= double(j+(Ny-2)*ColorModel.Dm->jproc())-double((Ny-2)*nprocy)*0.5;
double kglobal= double(k+(Nz-2)*ColorModel.Dm->kproc())-double((Nz-2)*nprocz)*0.5;
// Initialize phase position field for parallel bubble test
if ((iglobal*iglobal)+(jglobal*jglobal)+(kglobal*kglobal) < BubbleRadius*BubbleRadius){
ColorModel.Mask->id[n] = 2;
ColorModel.Mask->id[n] = 2;
}
else{
ColorModel.Mask->id[n]=1;
ColorModel.Mask->id[n]=1;
}
ColorModel.id[n] = ColorModel.Mask->id[n];
ColorModel.Dm->id[n] = ColorModel.Mask->id[n];
}
}
}
FILE *OUTFILE;
char LocalRankFilename[40];
sprintf(LocalRankFilename,"Bubble.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(ColorModel.id,1,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
// initialize the phase indicator field
}
int main( int argc, char **argv )
{
// Initialize
Utilities::startup( argc, argv );
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
// Initialize compute device
int device = ScaLBL_SetDevice( rank );
NULL_USE( device );
ScaLBL_DeviceBarrier();
comm.barrier();
Utilities::setErrorHandlers();
if ( rank == 0 ) {
printf( "********************************************************\n" );
printf( "Test Flow Adaptor with Color LBM \n" );
printf( "********************************************************\n" );
}
/* Populate the input database */
auto domain_db = std::make_shared<Database>();
auto color_db = std::make_shared<Database>();
auto vis_db = std::make_shared<Database>();
auto flow_db = std::make_shared<Database>();
auto analysis_db = std::make_shared<Database>();
auto db = std::make_shared<Database>();
domain_db->putVector<int>( "nproc", { 1, 1, 1 } );
domain_db->putVector<int>( "N", { 40, 40, 40 } );
domain_db->putVector<int>( "n", { 40, 40, 40 } );
domain_db->putScalar<int>( "BC", 0 );
flow_db->putScalar<double>("mass_fraction_factor",0.01);
flow_db->putScalar<int>("min_steady_timesteps",400);
flow_db->putScalar<int>("max_steady_timesteps",400);
flow_db->putScalar<int>("skiptimesteps",100);
color_db->putScalar<double>("alpha",0.01);
color_db->putScalar<int>("timestepMax",2000);
color_db->putVector<int>( "ComponentLabels", { 0, -1 } );
color_db->putVector<double>( "ComponentAffinity", { 0.0, 0.0 } );
db->putDatabase("Color",color_db);
db->putDatabase("Domain",domain_db);
db->putDatabase("FlowAdaptor",flow_db);
db->putDatabase("Visualization",vis_db);
db->putDatabase("Analysis",analysis_db);
ScaLBL_ColorModel ColorModel( rank, nprocs, comm );
ColorModel.color_db = color_db;
ColorModel.domain_db = domain_db;
ColorModel.vis_db = vis_db;
ColorModel.analysis_db = analysis_db;
ColorModel.db = db;
//ColorModel.ReadParams( filename );
ColorModel.SetDomain();
//ColorModel.ReadInput();
double radius=15.5;
InitializeBubble(ColorModel,radius);
ColorModel.Create(); // creating the model will create data structure to match the pore
// structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
double MLUPS=0.0;
int timestep = 0;
bool ContinueSimulation = true;
/* Variables for simulation protocols */
//auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "default" );
std::string PROTOCOL = "fractional flow";
/* 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 SEED_WATER = 0.01;
double ENDPOINT_THRESHOLD = 0.1;
double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled
MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 50000 );
ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
/* protocol specific key values */
if (PROTOCOL == "fractional flow")
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
else if (PROTOCOL == "seed water")
SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
/* analysis keys*/
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
}
/* Launch the simulation */
FlowAdaptor Adapt(ColorModel);
runAnalysis analysis(ColorModel);
while (ContinueSimulation){
/* this will run steady points */
if (PROTOCOL == "fractional flow" || PROTOCOL == "seed water" || PROTOCOL == "shell aggregation" || PROTOCOL == "image sequence" )
timestep += MAX_STEADY_TIME;
else
timestep += ColorModel.timestepMax;
/* Run the simulation timesteps*/
MLUPS = ColorModel.Run(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
if (ColorModel.timestep >= ColorModel.timestepMax){
ContinueSimulation = false;
}
/* Load a new image if image sequence is specified */
if (PROTOCOL == "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;
ContinueSimulation = false;
}
}
/*********************************************************/
/* update the fluid configuration with the flow adapter */
int skip_time = 0;
timestep = ColorModel.timestep;
/* get the averaged flow measures computed internally for the last simulation point*/
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;
double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M;
double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M;
double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M;
double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
/* stop simulation if previous point was sufficiently close to the endpoint*/
if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
if (ContinueSimulation && SKIP_TIMESTEPS > 0 ){
while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
timestep += ANALYSIS_INTERVAL;
if (PROTOCOL == "fractional flow") {
Adapt.UpdateFractionalFlow(ColorModel);
}
else if (PROTOCOL == "shell aggregation"){
double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange;
Adapt.ShellAggregation(ColorModel,target_volume_change);
}
else if (PROTOCOL == "seed water"){
Adapt.SeedPhaseField(ColorModel,SEED_WATER);
}
/* Run some LBM timesteps to let the system relax a bit */
MLUPS = ColorModel.Run(timestep);
/* Recompute the volume fraction now that the system has adjusted */
double volB = ColorModel.Averages->gwb.V;
double volA = ColorModel.Averages->gnb.V;
SaturationChange = volB/(volA + volB) - initialSaturation;
skip_time += ANALYSIS_INTERVAL;
}
if (rank==0) printf(" ********************************************************************* \n");
if (rank==0) printf(" Updated fractional flow with saturation change = %f \n", SaturationChange);
if (rank==0) printf(" Used protocol = %s \n", PROTOCOL.c_str());
if (rank==0) printf(" ********************************************************************* \n");
}
/*********************************************************/
}
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
return 0;
}

15
tests/lbpm_color_simulator.cpp
View File

@@ -92,6 +92,8 @@ int main( int argc, char **argv )
SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 50000 );
ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
/* protocol specific key values */
if (PROTOCOL == "image sequence" || PROTOCOL == "core flooding")
SKIP_TIMESTEPS = 0;
if (PROTOCOL == "fractional flow")
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
if (PROTOCOL == "seed water")
@@ -107,10 +109,14 @@ int main( int argc, char **argv )
runAnalysis analysis(ColorModel);
while (ContinueSimulation){
/* this will run steady points */
timestep += MAX_STEADY_TIME;
if (PROTOCOL == "fractional flow" || PROTOCOL == "seed water" || PROTOCOL == "shell aggregation" || PROTOCOL == "image sequence" )
timestep += MAX_STEADY_TIME;
else
timestep += ColorModel.timestepMax;
/* Run the simulation timesteps*/
MLUPS = ColorModel.Run(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
if (ColorModel.timestep > ColorModel.timestepMax){
if (ColorModel.timestep >= ColorModel.timestepMax){
ContinueSimulation = false;
}
@@ -175,10 +181,13 @@ int main( int argc, char **argv )
if (rank==0) printf(" ********************************************************************* \n");
}
/*********************************************************/
if (rank==0) printf(" (flatten density field) \n");
if (PROTOCOL == "fractional flow") {
Adapt.Flatten(ColorModel);
}
}
}
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );

1
threadpool/atomic_helpers.h
View File

@@ -222,7 +222,6 @@ inline int32_atomic atomic_fetch_and_or( int32_atomic volatile *v, int32_atomic
* \brief Fetch the current value and "ou" with given value
* \details Perform *v = (*v) | x, returning the previous value
* \return Returns the previous value before the "and" operation
* \param[in] v The pointer to the value to check and swap
* \param[in] v The pointer to the value to check and or
* \param[in] x The value to or
*/