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

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

Browse Source
This commit is contained in:
James E McClure 2018-07-30 15:50:13 -04:00
commit 5bd57c808b
11 changed files with 123 additions and 425 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
analysis/Minkowski.cpp
View File

@ -67,6 +67,22 @@ void Minkowski::Initialize()
vol_n = euler = Jn = An = Kn = 0.0;
}
double Minkowski::V(){
return vol_n_global;
}
double Minkowski::A(){
return An_global;
}
double Minkowski::J(){
return Jn_global;
}
double Minkowski::X(){
return euler_global;
}
void Minkowski::UpdateMeshValues()
{
@ -133,11 +149,9 @@ void Minkowski::ComputeLocal()
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm->id[n] != 0 ){
// 1-D index for this cube corner
if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
vol_n += 0.125;
}
// 1-D index for this cube corner
if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
vol_n += 0.125;
}
}

6
analysis/Minkowski.h
View File

@ -69,6 +69,11 @@ public:
DoubleArray SDn_y;
DoubleArray SDn_z;
double V();
double A();
double J();
double X();
//...........................................................................
Minkowski(std::shared_ptr <Domain> Dm);
~Minkowski();
@ -80,6 +85,7 @@ public:
void PrintAll();
int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
void SortBlobs();
};
#endif

2
common/ScaLBL.cpp
View File

@ -1121,14 +1121,12 @@ void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){
ScaLBL_DeviceBarrier();
//...................................................................................
// NOTE: AA Routine writes to opposite
// Unpack the gradributions on the device
//...................................................................................
//...Unpacking for x face(2,8,10,12,14)................................
ScaLBL_Gradient_Unpack(1.0,-1,0,0,dvcRecvDist_x,0,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,-1,0,dvcRecvDist_x,recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,1,0,dvcRecvDist_x,2*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,0,-1,dvcRecvDist_x,3*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,0,1,dvcRecvDist_x,4*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
//...................................................................................
//...Packing for X face(1,7,9,11,13)................................

14
example/Sph1896/input.db
View File

@ -1,3 +1,13 @@
MRT {
tau = 1.0
F = 0, 0, 1e-4
timestepMax = 1000
din = 1.0
dout = 1.0
Restart = false
flux = 0.0
}
Color {
tau = 1.0;
alpha = 1e-2;
@ -18,8 +28,8 @@ Color {
Domain {
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 80, 80, 80 // Size of local domain (Nx,Ny,Nz)
n_spheres = 1 // Number of spheres
n = 120, 120, 120 // Size of local domain (Nx,Ny,Nz)
nspheres = 1896 // Number of spheres
L = 1, 1, 1 // Length of domain (x,y,z)
BC = 0 // Boundary condition type
}

24
models/ColorModel.cpp
View File

@ -246,8 +246,8 @@ void ScaLBL_ColorModel::Create(){
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->last_interior; idx++){
if (idx == ScaLBL_Comm->next) idx = ScaLBL_Comm->first_interior;
for (int idx=0; idx<ScaLBL_Comm->LastInterior(); idx++){
if (idx == ScaLBL_Comm->LastExterior()) idx = ScaLBL_Comm->FirstInterior();
int n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n");
@ -316,8 +316,8 @@ void ScaLBL_ColorModel::Initialize(){
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->next, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
if (Dm->kproc()==0){
@ -365,9 +365,9 @@ void ScaLBL_ColorModel::Run(){
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@ -375,7 +375,7 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm_Regular->SendHalo(Phi);
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set BCs
@ -392,23 +392,23 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
@ -425,7 +425,7 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************

64
models/MRTModel.cpp
View File

@ -124,6 +124,7 @@ void ScaLBL_MRTModel::Create(){
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
}
void ScaLBL_MRTModel::Initialize(){
@ -162,21 +163,66 @@ void ScaLBL_MRTModel::Run(){
}
//************************************************************************/
stoptime = MPI_Wtime();
// cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl;
cputime = stoptime - starttime;
// cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 << " MLUPS" << endl;
double MLUPS = double(Np*timestep)/cputime/1000000;
// if (rank==0) printf("********************************************************\n");
// if (rank==0) printf("CPU time = %f \n", cputime);
// if (rank==0) printf("Lattice update rate (per process)= %f MLUPS \n", MLUPS);
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
}
void ScaLBL_MRTModel::VelocityField(double *Vz){
void ScaLBL_MRTModel::VelocityField(double *VELOCITY){
Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&Vz[0],&Velocity[2*Np],SIZE);
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
Morphology.Initialize();
Morphology.UpdateMeshValues();
Morphology.ComputeLocal();
Morphology.Reduce();
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
double mu = (tau-0.5)/3.f;
if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu,
Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
}

3
models/MRTModel.h
View File

@ -12,6 +12,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
class ScaLBL_MRTModel{
@ -54,6 +55,8 @@ public:
double *fq;
double *Velocity;
double *Pressure;
//Minkowski Morphology;
private:
MPI_Comm comm;

2
tests/CMakeLists.txt
View File

@ -2,7 +2,7 @@
#ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )

2
tests/GenerateSphereTest.cpp
View File

@ -494,7 +494,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
FILE *DIST = fopen(LocalRankFilename,"wb");
if (DIST==NULL) ERROR("Error opening file: ID.xxxxx");
fwrite(SignDist.data(),1,N,DIST);
fwrite(SignDist.data(),1,N*sizeof(double),DIST);
fclose(DIST);
//......................................................................

6
tests/TestPoiseuille.cpp
View File

@ -57,7 +57,7 @@ int main(int argc, char **argv)
MRT.Create(); // creating the model will create data structure to match the pore structure and allocate variables
MRT.Initialize(); // initializing the model will set initial conditions for variables
MRT.Run();
double *Vz; Vz= new double [MRT.Np];
double *Vz; Vz= new double [3*MRT.Np];
MRT.VelocityField(Vz);
if (rank == 0) printf("Force: %f,%f,%f \n",MRT.Fx,MRT.Fy,MRT.Fz);
@ -80,7 +80,7 @@ int main(int argc, char **argv)
n = MRT.Map(i,j,k);
//printf("%i,%i,%i; %i :",i,j,k,n);
if (n<0) {vz =0.f; printf(" b "); }
else { vz=Vz[n]; printf(" a "); }
else { vz=Vz[n+2*MRT.Np]; printf(" a "); }
printf("%f ",vz);
//Analytical solution
double x=1.f*i-1.5;
@ -97,7 +97,7 @@ int main(int argc, char **argv)
n = MRT.Map(i,j,k);
//printf("%i,%i,%i; %i :",i,j,k,n);
if (n<0) {vz =0.f; printf(" b "); }
else { vz=Vz[n]; printf(" a "); }
else { vz=Vz[n+2*MRT.Np]; printf(" a "); }
printf("%f ",vz);
//Analytical solution
double x=1.f*i-1.5;

401
tests/lbpm_permeability_simulator.cpp
View File

@ -10,7 +10,7 @@
#include "common/Communication.h"
#include "analysis/TwoPhase.h"
#include "common/MPI_Helpers.h"
#include "models/MRTModel.h"
//#define WRITE_SURFACES
/*
@ -42,396 +42,17 @@ int main(int argc, char **argv)
printf("Running Single Phase Permeability Calculation \n");
printf("********************************************************\n");
}
// Variables that specify the computational domain
string FILENAME;
int Nx,Ny,Nz; // local sub-domain size
int nspheres; // number of spheres in the packing
double Lx,Ly,Lz; // Domain length
double D = 1.0; // reference length for non-dimensionalization
// Color Model parameters
int timestepMax, interval;
double tau,Fx,Fy,Fz,tol,err;
double din,dout;
bool pBC,Restart;
int i,j,k,n;
int RESTART_INTERVAL=20000;
if (rank==0){
//.............................................................
// READ SIMULATION PARMAETERS FROM INPUT FILE
//.............................................................
ifstream input("Permeability.in");
// Line 1: model parameters (tau, alpha, beta, das, dbs)
input >> tau; // Viscosity parameter
// Line 2: External force components (Fx,Fy, Fz)
input >> Fx;
input >> Fy;
input >> Fz;
// Line 3: Pressure Boundary conditions
input >> Restart;
input >> pBC;
input >> din;
input >> dout;
// Line 4: time-stepping criteria
input >> timestepMax; // max no. of timesteps
input >> interval; // restart interval
input >> tol; // error tolerance
//.............................................................
//.......................................................................
// Reading the domain information file
//.......................................................................
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> Nx;
domain >> Ny;
domain >> Nz;
//domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
//.......................................................................
}
// **************************************************************
// Broadcast simulation parameters from rank 0 to all other procs
MPI_Barrier(comm);
//.................................................
MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
//MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm);
// MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm);
MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
MPI_Bcast(&interval,1,MPI_INT,0,comm);
MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm);
// Computational domain
MPI_Bcast(&Nx,1,MPI_INT,0,comm);
MPI_Bcast(&Ny,1,MPI_INT,0,comm);
MPI_Bcast(&Nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
//MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
MPI_Barrier(comm);
RESTART_INTERVAL=interval;
// **************************************************************
// **************************************************************
double rlx_setA = 1.f/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
if (nprocs != nprocx*nprocy*nprocz){
printf("nprocx = %i \n",nprocx);
printf("nprocy = %i \n",nprocy);
printf("nprocz = %i \n",nprocz);
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
if (rank==0){
printf("********************************************************\n");
printf("tau = %f \n", tau);
printf("Force(x) = %.5g \n", Fx);
printf("Force(y) = %.5g \n", Fy);
printf("Force(z) = %.5g \n", Fz);
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("********************************************************\n");
}
double viscosity=(tau-0.5)/3.0;
// Initialized domain and averaging framework for Two-Phase Flow
int BC=pBC;
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
Dm.CommInit();
TwoPhase Averages(Dm);
// Mask that excludes the solid phase
Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
MPI_Barrier(comm);
Nx += 2; Ny += 2; Nz += 2;
int N = Nx*Ny*Nz;
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
//.......................................................................
// Filenames used
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
char tmpstr[10];
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
// printf("Local File Name = %s \n",LocalRankFilename);
// .......... READ THE INPUT FILE .......................................
// char value;
char *id;
id = new char[N];
double sum, sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
//if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
double porosity, pore_vol;
//...........................................................................
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
//.......................................................................
// Read the signed distance
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
sum=0.f;
pore_vol = 0.0;
for ( k=0;k<Nz;k++){
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Averages.SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages.SDs(n) > 0.0){
sum++;
}
}
}
}
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("lbpm_permeability_simulator: Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_permeability_simulator: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
//.......................................................................
// Compute the media porosity, assign phase labels and solid composition
//.......................................................................
sum_local=0.0;
int Np=0; // number of local pore nodes
//.......................................................................
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (id[n] > 0){
sum_local+=1.0;
Np++;
}
}
}
}
MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm);
porosity = sum*iVol_global;
if (rank==0) printf("Media porosity = %f \n",porosity);
//.........................................................
// don't perform computations at the eight corners
id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
//.........................................................
MPI_Barrier(comm);
// Initialize communication structures in averaging domain
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
Mask.CommInit(comm);
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Mask);
// LBM variables
if (rank==0) printf ("Allocating distributions \n");
int Npad=(Np/16 + 2)*16;
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Mask.id,Np);
MPI_Barrier(comm);
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
int neighborSize=18*(Np*sizeof(int));
int *NeighborList;
// double *f_even,*f_odd;
double * dist;
double * Velocity;
double * Pressure;
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &dist, 19*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Pressure, 3*sizeof(double)*Np);
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
//...........................................................................
//...........................................................................
if (rank==0) printf("Setting the distributions, size = %i\n", N);
//...........................................................................
// Finalize setup for averaging domain
//Averages.SetupCubes(Dm);
Averages.UpdateSolid();
// Initialize two phase flow variables (all wetting phase)
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;
Averages.Phase(i,j,k) = -1.0;
Averages.SDn(i,j,k) = Averages.Phase(i,j,k);
Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k);
Averages.DelPhi(i,j,k) = 0.0;
Averages.Press(i,j,k) = 0.0;
Averages.Vel_x(i,j,k) = 0.0;
Averages.Vel_y(i,j,k) = 0.0;
Averages.Vel_z(i,j,k) = 0.0;
}
}
}
//.......................................................................
ScaLBL_D3Q19_Init(dist, Np);
int timestep = 0;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create and start timer............
double starttime,stoptime,cputime;
MPI_Barrier(comm);
starttime = MPI_Wtime();
//.........................................
double D32,Fo,Re,velocity,err1D,mag_force,vel_prev;
err = vel_prev = 1.0;
if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
//************ MAIN ITERATION LOOP ***************************************/
while (timestep < timestepMax && err > tol ){
timestep++;
ScaLBL_Comm.SendD3Q19AA(dist); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm.SendD3Q19AA(dist); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
if (timestep%500 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
ScaLBL_DeviceBarrier();
ScaLBL_D3Q19_Pressure(dist,Pressure,Np);
ScaLBL_D3Q19_Momentum(dist,Velocity,Np);
ScaLBL_Comm.RegularLayout(Map,Pressure,Averages.Press);
ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages.Vel_x);
ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages.Vel_y);
ScaLBL_Comm.RegularLayout(Map,&Velocity[2*Np],Averages.Vel_z);
// Way more work than necessary -- this is just to get the solid interfacial area!!
Averages.Initialize();
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
double vawx = Averages.vaw_global(0);
double vawy = Averages.vaw_global(1);
double vawz = Averages.vaw_global(2);
if (rank==0){
// ************* DIMENSIONLESS FORCHEIMER EQUATION *************************
// Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller
// Description of Non-Darcy Flows in Porous Medium Systems
// Physical Review E 87 (3), 033012
// Fo := density*D32^3*(density*force) / (viscosity^2)
// Re := density*D32*velocity / viscosity
// Fo = a*Re + b*Re^2
// *************************************************************************
//viscosity = (tau-0.5)*0.333333333333333333;
D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
printf("Sauter Mean Diameter = %f \n",D32);
mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
Fo = D32*D32*D32*mag_force/viscosity/viscosity;
// .... 1-D flow should be aligned with force ...
velocity = vawx*Fx/mag_force + vawy*Fy/mag_force + vawz*Fz/mag_force;
err1D = fabs(velocity-sqrt(vawx*vawx+vawy*vawy+vawz*vawz))/velocity;
//.......... Computation of the Reynolds number Re ..............
Re = D32*velocity/viscosity;
printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
printf("Relative error for 1D representation: %.5g \n",err1D);
printf("Dimensionless force: %5g \n", Fo);
printf("Reynolds number: %.5g \n", Re);
printf("Dimensionless Permeability (k/D^2): %.5g \n", Re/Fo);
}
}
}
//************************************************************************/
ScaLBL_DeviceBarrier();
MPI_Barrier(comm);
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
NULL_USE(RESTART_INTERVAL);
ScaLBL_MRTModel MRT(rank,nprocs,comm);
auto filename = argv[1];
MRT.ReadParams(filename);
MRT.SetDomain(); // this reads in the domain
MRT.ReadInput();
MRT.Create(); // creating the model will create data structure to match the pore structure and allocate variables
MRT.Initialize(); // initializing the model will set initial conditions for variables
MRT.Run();
double *Velocity; Velocity= new double [3*MRT.Np];
MRT.VelocityField(Velocity);
}
// ****************************************************
MPI_Barrier(comm);