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

Removed dead code from Domain.cpp, refactored lbpm_sphere_pp

Browse Source
This commit is contained in:
James E McClure
2018-05-16 09:21:06 -04:00
parent 9a22d71d0f
commit 17df5958c1
3 changed files with 258 additions and 474 deletions
Download Patch File Download Diff File Expand all files Collapse all files

410
common/Domain.cpp
View File

@@ -675,388 +675,6 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
}
/********************************************************
* Misc *
********************************************************/
double SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
/*
* This routine converts the data in the Distance array to a signed distance
* by solving the equation df/dt = sign(1-|grad f|), where Distance provides
* the values of f on the mesh associated with domain Dm
* It has been tested with segmented data initialized to values [-1,1]
* and will converge toward the signed distance to the surface bounding the associated phases
*/
int Q=26;
int q,i,j,k;
double dt=0.1;
int in,jn,kn;
double Dqx,Dqy,Dqz,Dx,Dy,Dz,W;
double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm;
double TotalVariation=0.0;
const static int D3Q27[26][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
{1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1},{1,1,1},{-1,-1,-1},{1,1,-1},{-1,-1,1},
{-1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1}};
double weights[26];
// Compute the weights from the finite differences
for (q=0; q<Q; q++){
weights[q] = sqrt(1.0*(D3Q27[q][0]*D3Q27[q][0]) + 1.0*(D3Q27[q][1]*D3Q27[q][1]) + 1.0*(D3Q27[q][2]*D3Q27[q][2]));
}
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
int count = 0;
while (count < timesteps){
// Communicate the halo of values
fillData.fill(Distance);
TotalVariation=0.0;
// Execute the next timestep
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
//sign = Distance(i,j,k) / fabs(Distance(i,j,k));
sign = -1;
if (ID[n] == 1) sign = 1;
/*
if (!(i+1<Nx)) nx=0.5*Distance(i,j,k);
else nx=0.5*Distance(i+1,j,k);;
if (!(j+1<Ny)) ny=0.5*Distance(i,j,k);
else ny=0.5*Distance(i,j+1,k);
if (!(k+1<Nz)) nz=0.5*Distance(i,j,k);
else nz=0.5*Distance(i,j,k+1);
if (i<1) nx-=0.5*Distance(i,j,k);
else nx-=0.5*Distance(i-1,j,k);
if (j<1) ny-=0.5*Distance(i,j,k);
else ny-=0.5*Distance(i,j-1,k);
if (k<1) nz-=0.5*Distance(i,j,k);
else nz-=0.5*Distance(i,j,k-1);
*/
//............Compute the Gradient...................................
nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
W = 0.0; Dx = Dy = Dz = 0.0;
// also ignore places where the gradient is zero since this will not
// result in any local change to Distance
if (nx*nx+ny*ny+nz*nz > 0.0 ){
for (q=0; q<26; q++){
Cqx = 1.0*D3Q27[q][0];
Cqy = 1.0*D3Q27[q][1];
Cqz = 1.0*D3Q27[q][2];
// get the associated neighbor
in = i + D3Q27[q][0];
jn = j + D3Q27[q][1];
kn = k + D3Q27[q][2];
// make sure the neighbor is in the domain (periodic BC)
/* if (in < 0 ) in +=Nx;
* don't need this in parallel since MPI handles the halos
if (jn < 0 ) jn +=Ny;
if (kn < 0 ) kn +=Nz;
if (!(in < Nx) ) in -=Nx;
if (!(jn < Ny) ) jn -=Ny;
if (!(kn < Nz) ) kn -=Nz;
// symmetric boundary
if (in < 0 ) in = i;
if (jn < 0 ) jn = j;
if (kn < 0 ) kn = k;
if (!(in < Nx) ) in = i;
if (!(jn < Ny) ) jn = k;
if (!(kn < Nz) ) kn = k;
*/
// Compute the gradient using upwind finite differences
Dqx = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqx;
Dqy = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqy;
Dqz = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqz;
// Only include upwind derivatives
if (sign*(nx*Cqx + ny*Cqy + nz*Cqz) < 0.0 ){
Dx += Dqx;
Dy += Dqy;
Dz += Dqz;
W += weights[q];
}
}
// Normalize by the weight to get the approximation to the gradient
if (fabs(W) > 0.0){
Dx /= W;
Dy /= W;
Dz /= W;
}
norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
}
else{
norm = 0.0;
}
Distance(i,j,k) += dt*sign*(1.0 - norm);
TotalVariation += dt*sign*(1.0 - norm);
// Disallow any change in phase
// if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k);
}
}
}
TotalVariation /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2);
count++;
}
return TotalVariation;
}
void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad)
{
// Read in the full sphere pack
//...... READ IN THE SPHERES...................................
cout << "Reading the packing file..." << endl;
FILE *fid = fopen("pack.out","rb");
INSIST(fid!=NULL,"Error opening pack.out");
//.........Trash the header lines..........
char line[100];
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
//........read the spheres..................
// We will read until a blank like or end-of-file is reached
int count = 0;
while ( !feof(fid) && fgets(line,100,fid)!=NULL ) {
char* line2 = line;
List_cx[count] = strtod(line2,&line2);
List_cy[count] = strtod(line2,&line2);
List_cz[count] = strtod(line2,&line2);
List_rad[count] = strtod(line2,&line2);
count++;
}
cout << "Number of spheres extracted is: " << count << endl;
INSIST( count==nspheres, "Specified number of spheres is probably incorrect!" );
// .............................................................
}
void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
char value;
int N = Nx*Ny*Nz; // Domain size, including the halo
double hx,hy,hz;
double x,y,z;
double cx,cy,cz,r;
int imin,imax,jmin,jmax,kmin,kmax;
int p,i,j,k,n;
//............................................
double min_x,min_y,min_z;
// double max_x,max_y,max_z;
//............................................
// Lattice spacing for the entire domain
// It should generally be true that hx=hy=hz
// Otherwise, you will end up with ellipsoids
hx = Lx/(Nx*nprocx-1);
hy = Ly/(Ny*nprocy-1);
hz = Lz/(Nz*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*Nx-1)*hx;
min_y = double(jproc*Ny-1)*hy;
min_z = double(kproc*Nz-1)*hz;
// max_x = ((iproc+1)*Nx+1)*hx;
// max_y = ((jproc+1)*Ny+1)*hy;
// max_z = ((kproc+1)*Nz+1)*hz;
//............................................
//............................................
// Pre-initialize local ID
for (n=0;n<N;n++){
ID[n]=1;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-r)/hx)-1;
imax = int ((cx+r)/hx)+1;
jmin = int ((cy-r)/hy)-1;
jmax = int ((cy+r)/hy)+1;
kmin = int ((cz-r)/hz)-1;
kmax = int ((cz+r)/hz)+1;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) imin = Nx;
if (imax<0) imax = 0;
if (imax>Nx) imax = Nx;
if (jmin<0) jmin = 0;
if (jmin>Ny) jmin = Ny;
if (jmax<0) jmax = 0;
if (jmax>Ny) jmax = Ny;
if (kmin<0) kmin = 0;
if (kmin>Nz) kmin = Nz;
if (kmax<0) kmax = 0;
if (kmax>Nz) kmax = Nz;
// Loop over the domain for this sphere (may be null)
for (i=imin;i<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// Initialize ID value to 'fluid (=1)'
x = i*hx;
y = j*hy;
z = k*hz;
value = 1;
// if inside sphere, set to zero
if ( (cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z) < r*r){
value=0;
}
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
if ( ID[n] != 0 ){
ID[n] = value;
}
}
}
}
}
}
void SignedDistance(double *Distance, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
int N = Nx*Ny*Nz; // Domain size, including the halo
double hx,hy,hz;
double x,y,z;
double cx,cy,cz,r;
int imin,imax,jmin,jmax,kmin,kmax;
int p,i,j,k,n;
//............................................
double min_x,min_y,min_z;
double distance;
//............................................
// Lattice spacing for the entire domain
// It should generally be true that hx=hy=hz
// Otherwise, you will end up with ellipsoids
hx = Lx/((Nx-2)*nprocx-1);
hy = Ly/((Ny-2)*nprocy-1);
hz = Lz/((Nz-2)*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*(Nx-2)-1)*hx;
min_y = double(jproc*(Ny-2)-1)*hy;
min_z = double(kproc*(Nz-2)-1)*hz;
//............................................
//............................................
// Pre-initialize Distance
for (n=0;n<N;n++){
Distance[n]=100.0;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-2*r)/hx);
imax = int ((cx+2*r)/hx)+2;
jmin = int ((cy-2*r)/hy);
jmax = int ((cy+2*r)/hy)+2;
kmin = int ((cz-2*r)/hz);
kmax = int ((cz+2*r)/hz)+2;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) imin = Nx;
if (imax<0) imax = 0;
if (imax>Nx) imax = Nx;
if (jmin<0) jmin = 0;
if (jmin>Ny) jmin = Ny;
if (jmax<0) jmax = 0;
if (jmax>Ny) jmax = Ny;
if (kmin<0) kmin = 0;
if (kmin>Nz) kmin = Nz;
if (kmax<0) kmax = 0;
if (kmax>Nz) kmax = Nz;
// Loop over the domain for this sphere (may be null)
for (i=imin;i<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// x,y,z is distance in physical units
x = i*hx;
y = j*hy;
z = k*hz;
// if inside sphere, set to zero
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
// Compute the distance
distance = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z)) - r;
// Assign the minimum distance
if (distance < Distance[n]) Distance[n] = distance;
}
}
}
}
// Map the distance to lattice units
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
}
void WriteLocalSolidID(char *FILENAME, char *ID, int N)
{
char value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = ID[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N)
{
double value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = Distance[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np)
{
int q,n;
@@ -1078,18 +696,15 @@ void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq
}
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np)
void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np)
{
int q=0, n=0;
double value=0;
ifstream File(FILENAME,ios::binary);
for (n=0; n<Np; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
cDen[n] = value;
File.read((char*) &value, sizeof(value));
cDen[Np+n] = value;
// Read the even distributions
cPhi[n] = value;
// Read the distributions
for (q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value;
@@ -1098,24 +713,5 @@ void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np)
File.close();
}
void ReadBinaryFile(char *FILENAME, double *Data, int N)
{
int n;
double value;
ifstream File(FILENAME,ios::binary);
if (File.good()){
for (n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Data[n] = value;
}
}
else {
for (n=0; n<N; n++) Data[n] = 1.2e-34;
}
File.close();
}

17
common/Domain.h
View File

@@ -230,26 +230,9 @@ private:
TYPE *d_gcw;
};
void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad);
void SignedDistance(double *Distance, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz);
void WriteLocalSolidID(char *FILENAME, char *ID, int N);
void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N);
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np);
void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np);
void ReadBinaryFile(char *FILENAME, double *Data, int N);
#endif

305
tests/lbpm_sphere_pp.cpp
View File

@@ -19,6 +19,244 @@
using namespace std;
void WriteLocalSolidID(char *FILENAME, char *ID, int N)
{
char value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = ID[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void WriteLocalSolidDistance(char *FILENAME, double *Distance, int N)
{
double value;
ofstream File(FILENAME,ios::binary);
for (int n=0; n<N; n++){
value = Distance[n];
File.write((char*) &value, sizeof(value));
}
File.close();
}
void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad)
{
// Read in the full sphere pack
//...... READ IN THE SPHERES...................................
cout << "Reading the packing file..." << endl;
FILE *fid = fopen("pack.out","rb");
INSIST(fid!=NULL,"Error opening pack.out");
//.........Trash the header lines..........
char line[100];
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
fgetl(line, 100, fid);
//........read the spheres..................
// We will read until a blank like or end-of-file is reached
int count = 0;
while ( !feof(fid) && fgets(line,100,fid)!=NULL ) {
char* line2 = line;
List_cx[count] = strtod(line2,&line2);
List_cy[count] = strtod(line2,&line2);
List_cz[count] = strtod(line2,&line2);
List_rad[count] = strtod(line2,&line2);
count++;
}
cout << "Number of spheres extracted is: " << count << endl;
INSIST( count==nspheres, "Specified number of spheres is probably incorrect!" );
// .............................................................
}
void AssignLocalSolidID(char *ID, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
char value;
int N = Nx*Ny*Nz; // Domain size, including the halo
double hx,hy,hz;
double x,y,z;
double cx,cy,cz,r;
int imin,imax,jmin,jmax,kmin,kmax;
int p,i,j,k,n;
//............................................
double min_x,min_y,min_z;
// double max_x,max_y,max_z;
//............................................
// Lattice spacing for the entire domain
// It should generally be true that hx=hy=hz
// Otherwise, you will end up with ellipsoids
hx = Lx/(Nx*nprocx-1);
hy = Ly/(Ny*nprocy-1);
hz = Lz/(Nz*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*Nx-1)*hx;
min_y = double(jproc*Ny-1)*hy;
min_z = double(kproc*Nz-1)*hz;
// max_x = ((iproc+1)*Nx+1)*hx;
// max_y = ((jproc+1)*Ny+1)*hy;
// max_z = ((kproc+1)*Nz+1)*hz;
//............................................
//............................................
// Pre-initialize local ID
for (n=0;n<N;n++){
ID[n]=1;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-r)/hx)-1;
imax = int ((cx+r)/hx)+1;
jmin = int ((cy-r)/hy)-1;
jmax = int ((cy+r)/hy)+1;
kmin = int ((cz-r)/hz)-1;
kmax = int ((cz+r)/hz)+1;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) imin = Nx;
if (imax<0) imax = 0;
if (imax>Nx) imax = Nx;
if (jmin<0) jmin = 0;
if (jmin>Ny) jmin = Ny;
if (jmax<0) jmax = 0;
if (jmax>Ny) jmax = Ny;
if (kmin<0) kmin = 0;
if (kmin>Nz) kmin = Nz;
if (kmax<0) kmax = 0;
if (kmax>Nz) kmax = Nz;
// Loop over the domain for this sphere (may be null)
for (i=imin;i<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// Initialize ID value to 'fluid (=1)'
x = i*hx;
y = j*hy;
z = k*hz;
value = 1;
// if inside sphere, set to zero
if ( (cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z) < r*r){
value=0;
}
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
if ( ID[n] != 0 ){
ID[n] = value;
}
}
}
}
}
}
void SignedDistance(double *Distance, int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad,
double Lx, double Ly, double Lz, int Nx, int Ny, int Nz,
int iproc, int jproc, int kproc, int nprocx, int nprocy, int nprocz)
{
// Use sphere lists to determine which nodes are in porespace
// Write out binary file for nodes
int N = Nx*Ny*Nz; // Domain size, including the halo
double hx,hy,hz;
double x,y,z;
double cx,cy,cz,r;
int imin,imax,jmin,jmax,kmin,kmax;
int p,i,j,k,n;
//............................................
double min_x,min_y,min_z;
double distance;
//............................................
// Lattice spacing for the entire domain
// It should generally be true that hx=hy=hz
// Otherwise, you will end up with ellipsoids
hx = Lx/((Nx-2)*nprocx-1);
hy = Ly/((Ny-2)*nprocy-1);
hz = Lz/((Nz-2)*nprocz-1);
//............................................
// Get maximum and minimum for this domain
// Halo is included !
min_x = double(iproc*(Nx-2)-1)*hx;
min_y = double(jproc*(Ny-2)-1)*hy;
min_z = double(kproc*(Nz-2)-1)*hz;
//............................................
//............................................
// Pre-initialize Distance
for (n=0;n<N;n++){
Distance[n]=100.0;
}
//............................................
//............................................
// .........Loop over the spheres.............
for (p=0;p<nspheres;p++){
// Get the sphere from the list, map to local min
cx = List_cx[p] - min_x;
cy = List_cy[p] - min_y;
cz = List_cz[p] - min_z;
r = List_rad[p];
// Check if
// Range for this sphere in global indexing
imin = int ((cx-2*r)/hx);
imax = int ((cx+2*r)/hx)+2;
jmin = int ((cy-2*r)/hy);
jmax = int ((cy+2*r)/hy)+2;
kmin = int ((cz-2*r)/hz);
kmax = int ((cz+2*r)/hz)+2;
// Obviously we have to do something at the edges
if (imin<0) imin = 0;
if (imin>Nx) imin = Nx;
if (imax<0) imax = 0;
if (imax>Nx) imax = Nx;
if (jmin<0) jmin = 0;
if (jmin>Ny) jmin = Ny;
if (jmax<0) jmax = 0;
if (jmax>Ny) jmax = Ny;
if (kmin<0) kmin = 0;
if (kmin>Nz) kmin = Nz;
if (kmax<0) kmax = 0;
if (kmax>Nz) kmax = Nz;
// Loop over the domain for this sphere (may be null)
for (i=imin;i<imax;i++){
for (j=jmin;j<jmax;j++){
for (k=kmin;k<kmax;k++){
// x,y,z is distance in physical units
x = i*hx;
y = j*hy;
z = k*hz;
// if inside sphere, set to zero
// get the position in the list
n = k*Nx*Ny+j*Nx+i;
// Compute the distance
distance = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z)) - r;
// Assign the minimum distance
if (distance < Distance[n]) Distance[n] = distance;
}
}
}
}
// Map the distance to lattice units
for (n=0; n<N; n++) Distance[n] = Distance[n]/hx;
}
int main(int argc, char **argv)
{
@@ -52,53 +290,27 @@ int main(int argc, char **argv)
printf("********************************************************\n");
}
// Variables that specify the computational domain
string FILENAME;
unsigned int nBlocks, nthreads;
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
// Load inputs
string FILENAME = argv[1];
// Load inputs
if (rank==0) printf("Loading input database \n");
auto db = std::make_shared<Database>( FILENAME );
auto domain_db = db->getDatabase( "Domain" );
int Nx = domain_db->getVector<int>( "n" )[0];
int Ny = domain_db->getVector<int>( "n" )[1];
int Nz = domain_db->getVector<int>( "n" )[2];
int nprocx = domain_db->getVector<int>( "nproc" )[0];
int nprocy = domain_db->getVector<int>( "nproc" )[1];
int nprocz = domain_db->getVector<int>( "nproc" )[2];
int nspheres = domain_db->getScalar<int>( "nsphere" );
int Lx = domain_db->getVector<double>( "L" )[0];
int Ly = domain_db->getVector<double>( "L" )[1];
int Lz = domain_db->getVector<double>( "L" )[2];
int i,j,k,n;
if (rank==0){
//.......................................................................
// 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);
//.................................................
// 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);
// **************************************************************
if (nprocs != nprocx*nprocy*nprocz){
printf("nprocx = %i \n",nprocx);
printf("nprocy = %i \n",nprocy);
@@ -106,13 +318,6 @@ int main(int argc, char **argv)
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz,
rank_yz, rank_YZ, rank_yZ, rank_Yz );
MPI_Barrier(comm);
Nz += 2;
Nx = Ny = Nz; // Cubic domain