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:
James E McClure 2016-11-24 10:51:13 -05:00
commit ac1222ee81
7 changed files with 945 additions and 455 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
tests/CMakeLists.txt
View File

@ -1,5 +1,6 @@
# Copy files for the tests
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )

26
tests/lbpm_morph_pp.cpp
View File

@ -233,30 +233,36 @@ int main(int argc, char **argv)
if (rank==0){
PORESIZE=fopen("PoreSize.hist","w");
printf(" writing PoreSize.hist \n");
int PoreCount=0;
int Count;
uint64_t PoreCount=0;
uint64_t Count;
double PoreVol=0.f;
for (int idx=0; idx<NumBins; idx++){
double BinCenter=MinPoreSize+idx*BinWidth;
Count=GlobalHistogram[idx];
PoreCount+=Count;
fprintf(PORESIZE,"%i %f\n",Count,BinCenter);
PoreVol+=Count*BinCenter*BinCenter*BinCenter;
fprintf(PORESIZE,"%lu %f\n",Count,BinCenter);
}
fclose(PORESIZE);
// Compute quartiles
//printf("Total pores: %lu\n",PoreCount);
double Q1,Q2,Q3,Q4;
int Qval=PoreCount/4;
//uint64_t Qval=PoreCount/4;
double Qval=PoreVol*0.25;
Q1=Q2=Q3=MinPoreSize;
Q4=MaxPoreSize;
Count=0;
//printf("Volume per quartile %f\n",Qval);
PoreVol=0.f;
for (int idx=0; idx<NumBins; idx++){
double BinCenter=MinPoreSize+idx*BinWidth;
Count+=GlobalHistogram[idx];
if (Count<Qval) Q1+=BinWidth;
if (Count<2*Qval) Q2+=BinWidth;
if (Count<3*Qval) Q3+=BinWidth;
Count=GlobalHistogram[idx];
PoreVol+=Count*BinCenter*BinCenter*BinCenter;
if (PoreVol<Qval) Q1+=BinWidth;
if (PoreVol<2*Qval) Q2+=BinWidth;
if (PoreVol<3*Qval) Q3+=BinWidth;
}
printf("Quartiles for pore size distribution \n");
printf("Quartiles (volumetric) for pore size distribution \n");
printf("Q1 %f\n",Q1);
printf("Q2 %f\n",Q2);
printf("Q3 %f\n",Q3);

123
tests/lbpm_morphdrain_pp.cpp
View File

@ -167,43 +167,6 @@ int main(int argc, char **argv)
fclose(IDFILE);
int count,countGlobal,totalGlobal;
count = 0;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
else{
// initially saturated with wetting phase
//id[n] = 2;
count++;
}
}
}
}
// total Global is the number of nodes in the pore-space
MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float porosity=float(totalGlobal)/(nprocx*nprocy*nprocz*(nx-2)*(ny-2)*(nz-2));
if (rank==0) printf("Media Porosity: %f \n",porosity);
count = 0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
n=k*Nx*Ny+j*Nx+i;
if (id[n] == 2){
count++;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sw= double(countGlobal)/totalGlobal;
if (rank==0) printf("Initial saturation (from ID.xxxxx files)=%f\n",sw)
Dm.CommInit(comm);
int iproc = Dm.iproc;
int jproc = Dm.jproc;
@ -269,10 +232,30 @@ int main(int argc, char **argv)
int Nx = nx;
int Ny = ny;
int Nz = nz;
double sw=1.f;
int GlobalNumber = 1;
double radius,Rcrit;
int count,countGlobal,totalGlobal;
count = 0;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
else{
// initially saturated with wetting phase
id[n] = 2;
count++;
}
}
}
}
// total Global is the number of nodes in the pore-space
MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
float porosity=float(totalGlobal)/(nprocx*nprocy*nprocz*(nx-2)*(ny-2)*(nz-2));
if (rank==0) printf("Media Porosity: %f \n",porosity);
double radius,Rcrit_new;
radius = 0.0;
// Layer the inlet with NWP
if (kproc == 0){
@ -281,27 +264,35 @@ int main(int argc, char **argv)
n = j*nx+i;
// n = nx*ny + j*nx+i;
id[n]=1;
if (SignDist(i,j,k) > radius){
radius=SignDist(i,j,k);
if (SignDist(i,j,0) > radius){
radius=SignDist(i,j,0);
}
}
}
}
MPI_Allreduce(&radius,&Rcrit,1,MPI_DOUBLE,MPI_MAX,comm);
int Window=int(Rcrit);
MPI_Allreduce(&radius,&Rcrit_new,1,MPI_DOUBLE,MPI_MAX,comm);
if (rank==0) printf("Starting morhpological drainage with critical radius = %f \n",Rcrit);
if (rank==0) printf("Starting morhpological drainage with critical radius = %f \n",Rcrit_new);
int imin,jmin,kmin,imax,jmax,kmax;
// Decrease the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
while (sw > SW){
double deltaR=0.01; // amount to change the radius in voxel units
double Rcrit_old;
double sw_old=1.0; // initial WP saturation set to one
double sw_new=1.0; // initial WP saturation set to one
double sw_diff_old = 1.0;
double sw_diff_new = 1.0;
while (sw_new > SW){
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR;// decrease critical radius
sw_old = sw_new;
sw_diff_old = sw_diff_new;
// decrease critical radius
Rcrit -= deltaR;
Window=int(Rcrit);
int Window=round(Rcrit_new);
GlobalNumber = 1;
while (GlobalNumber != 0){
@ -312,7 +303,7 @@ int main(int argc, char **argv)
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
n = k*nx*ny + j*nx+i;
if (id[n] == 1 && SignDist(i,j,k) > Rcrit){
if (id[n] == 1 && SignDist(i,j,k) > Rcrit_new){
// loop over the window and update
imin=max(1,i-Window);
jmin=max(1,j-Window);
@ -325,7 +316,7 @@ int main(int argc, char **argv)
for (ii=imin; ii<imax; ii++){
int nn = kk*nx*ny+jj*nx+ii;
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (id[nn] == 2 && dsq <= Rcrit*Rcrit){
if (id[nn] == 2 && dsq <= Rcrit_new*Rcrit_new){
LocalNumber++;
id[nn]=1;
}
@ -454,10 +445,36 @@ int main(int argc, char **argv)
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sw= double(countGlobal)/totalGlobal;
if (rank==0) printf("Final saturation=%f\n",sw);
sw_new= double(countGlobal)/totalGlobal;
sw_diff_new = abs(sw_new-SW);
// if (rank==0){
// printf("Final saturation=%f\n",sw_new);
// printf("Final critical radius=%f\n",Rcrit_new);
//
// }
}
if (sw_diff_new<sw_diff_old){
if (rank==0){
printf("Final saturation=%f\n",sw_new);
printf("Final critical radius=%f\n",Rcrit_new);
}
}
else{
if (rank==0){
printf("Final saturation=%f\n",sw_old);
printf("Final critical radius=%f\n",Rcrit_old);
}
}
// if (rank==0){
// printf("Final saturation=%f\n",sw);
// printf("Final critical radius=%f\n",Rcrit);
//
// }
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);

85
tests/lbpm_morphopen_pp.cpp
View File

@ -54,7 +54,8 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
double Rcrit=0.f;
//double Rcrit_new=1.f; // Hard-coded 'Rcrit' to avoid any calculations under resolutions.
double Rcrit_new=0.f;
double SW=strtod(argv[1],NULL);
if (rank==0){
//printf("Critical radius %f (voxels)\n",Rcrit);
@ -330,39 +331,35 @@ int main(int argc, char **argv)
int Nx = nx;
int Ny = ny;
int Nz = nz;
double sw = 0.f;
int GlobalNumber = 1;
double sw_old=1.0;
double sw_new=1.0;
double sw_diff_old = 1.0;
double sw_diff_new = 1.0;
int imin,jmin,kmin,imax,jmax,kmax;
// Increase the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
while (sw<SW)
{
double Rcrit_old;
Rcrit += deltaR;
int Window=round(Rcrit);
int GlobalNumber = 1;
int imin,jmin,kmin,imax,jmax,kmax;
Rcrit_new = maxdistGlobal;
while (sw_new > SW)
{
sw_diff_old = sw_diff_new;
sw_old = sw_new;
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR;
int Window=round(Rcrit_new);
if (Window == 0) Window = 1; // If Window = 0 at the begining, after the following process will have sw=1.0
// and sw<Sw will be immediately broken
int LocalNumber=0;
// Initialization: saturate medium with wetting phase - need this for each iteraction before SW is met
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
else{
// initially saturated with wetting phase
id[n] = 2;
}
}
}
}
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;
if (SignDist(i,j,k) > Rcrit){
if (SignDist(i,j,k) > Rcrit_new){
// loop over the window and update
imin=max(1,i-Window);
jmin=max(1,j-Window);
@ -375,13 +372,14 @@ int main(int argc, char **argv)
for (ii=imin; ii<imax; ii++){
int nn = kk*nx*ny+jj*nx+ii;
double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
if (id[nn] == 2 && dsq <= Rcrit*Rcrit){
if (id[nn] == 2 && dsq <= Rcrit_new*Rcrit_new){
LocalNumber++;
id[nn]=1;
}
}
}
}
}
// move on
}
@ -480,23 +478,34 @@ int main(int argc, char **argv)
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sw = float(countGlobal)/totalGlobal;
if (rank==0)
{
printf("Final saturation=%f\n",sw);
printf("Final critical radius=%f\n",Rcrit);
sw_new = float(countGlobal)/totalGlobal;
sw_diff_new = abs(sw_new-SW);
// for test only
// if (rank==0){
// printf("Final saturation=%f\n",sw_new);
// printf("Final critical radius=%f\n",Rcrit_new);
//
// }
}
if (sw_diff_new<sw_diff_old){
if (rank==0){
printf("Final saturation=%f\n",sw_new);
printf("Final critical radius=%f\n",Rcrit_new);
}
}
// Restore the solid phase
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
n = k*nx*ny+j*nx+i;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
}
}
}
else{
if (rank==0){
printf("Final saturation=%f\n",sw_old);
printf("Final critical radius=%f\n",Rcrit_old);
}
}
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);

591
tests/lbpm_nondarcy_simulator.cpp Normal file
View File

@ -0,0 +1,591 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/TwoPhase.h"
#include "common/MPI_Helpers.h"
//#define WRITE_SURFACES
/*
* Simulator for two-phase flow in porous media
* James E. McClure 2013-2014
*/
using namespace std;
//*************************************************************************
// Steady State Single-Phase LBM to generate non-Darcy curves
//*************************************************************************
int main(int argc, char **argv)
{
std::string help("--help");
std::string arg1("");
if (argc > 1) arg1=argv[1];
if (help.compare(arg1) == 0){
printf("********************************************************** \n");
printf("Pore-scale lattice Boltzmann simulator for non-Darcy flow in porous media \n \n");
printf("Simulate non-Darcy flow in porous media \n");
printf(" MPI-based lattice Boltzmann simulator \n");
printf(" Multi-relaxation time (MRT) D3Q19 \n \n");
printf(" Launch with MPI (e.g.) \n \n");
printf(" mpirun -np $NUMPROCS $LBPM_WIA_DIR/lbpm_nondarcy_simulator \n \n");
printf("**********CITATION********** \n");
printf(" Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller\n");
printf(" Description of Non-Darcy Flows in Porous Medium Systems \n");
printf(" Physical Review E 87 (3), 033012 \n \n");
printf("**********INPUT********** \n");
printf("1. Domain.in (describes the simulation domain and domain decomposition) \n");
printf(" ----(e.g. Domain.in)-----\n");
printf(" nprocx nprocy nprocz (process grid)\n");
printf(" Nx Ny Nz (local sub-domain)\n");
printf(" Lx Ly Lz (physical domain size) \n");
printf(" --------------------------\n");
printf("2. SignDist.xxxxx (Distance map of porespace) \n");
printf(" - one file for each MPI process \n");
printf(" - double precision values \n");
printf(" - dimensions are [Nx+2,Ny+2,Nz+2] (include halo)\n \n");
//printf("3. parameters for LBM are hard-coded! \n \n");
printf("**********OUTPUT********** \n");
printf("1. nondary.csv - list of averaged quantities obtained from steady-state flow fields\n");
printf(" - D32 - Sauter mean grain diamter \n");
printf(" - vx - average velocity in the x-direction \n");
printf(" - vy - average velocity in the y-direction \n");
printf(" - vz - average velocity in the z-direction \n");
printf(" - Fx - body force applied in the x-direction \n");
printf(" - Fy - body force applied in the y-direction \n");
printf(" - Fz - body force applied in the z-direction \n");
printf(" - Fo - Gallilei number \n");
printf(" - Re - Reynolds number \n");
printf("********************************************************** \n");
/*printf("*******DIMENSIONLESS FORCHEIMER EQUATION********\n");
printf(" - force = |F| (total magnitude of force) \n");
printf(" - velocity = F dot v / |F| (flow velocity aligned with force) \n");
printf(" - Fo = density*D32^3*(density*force) / (viscosity^2) \n");
printf(" - Re = density*D32*velocity / viscosity \n");
printf(" - Fo = a*Re + b*Re^2 \n");
*/
// *************************************************************************
}
else {
//*****************************************
// ***** MPI STUFF ****************
//*****************************************
// Initialize MPI
int rank,nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
//*****************************************
// MPI ranks for all 18 neighbors
//**********************************
int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
int rank_xy,rank_XY,rank_xY,rank_Xy;
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//**********************************
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
double REYNOLDS_NUMBER = 100.f;
if (argc > 1){
REYNOLDS_NUMBER=strtod(argv[1],NULL);
}
if (rank == 0){
printf("********************************************************\n");
printf("Simulating Single Phase Non-Darcy Curve, Re < %f \n",REYNOLDS_NUMBER);
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){
//.......................................................................
// 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;
//.......................................................................
/*
* Set simulation parameters internally
*/
tau=1.f;
Fx = 0.f;
Fy = 0.f;
Fz = 1.0e-7;
pBC = 0;
din = 1.0;
dout = 1.0;
timestepMax = nprocz*Nz*100;
interval = 500;
tol = 1.0e-4;
}
// **************************************************************
// 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 rlxA = 1.f/tau;
double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
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) = %f \n", Fx);
printf("Force(y) = %f \n", Fy);
printf("Force(z) = %f \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);
TwoPhase Averages(Dm);
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);
Nx += 2; Ny += 2; Nz += 2;
int N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
//.......................................................................
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];
int sum = 0;
double sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (pBC) 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;
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
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++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 0;
}
}
}
sum=0;
pore_vol = 0.0;
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 (Averages.SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages.SDs(n) > 0.0){
sum++;
}
}
}
}
// Set up kstart, kfinish so that the reservoirs are excluded from averaging
int kstart,kfinish;
kstart = 1;
kfinish = Nz-1;
if (pBC && kproc==0) kstart = 4;
if (pBC && kproc==nprocz-1) kfinish = Nz-4;
// Compute the pore volume
sum_local = 0.0;
for ( k=kstart;k<kfinish;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;
}
}
}
}
MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
// MPI_Allreduce(&sum_local,&porosity,1,MPI_DOUBLE,MPI_SUM,comm);
porosity = pore_vol*iVol_global;
if (rank==0) printf("Media porosity = %f \n",porosity);
//.........................................................
// If pressure boundary conditions are applied remove solid
if (pBC && kproc == 0){
for (k=0; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
}
}
}
}
if (pBC && kproc == nprocz-1){
for (k=Nz-3; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
id[n] = 2;
Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
}
}
}
}
//.........................................................
// 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;
//.........................................................
// Initialize communication structures in averaging domain
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
Dm.CommInit(comm);
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Dm);
//...........device phase ID.................................................
if (rank==0) printf ("Copying phase ID to device \n");
char *ID;
AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
// Copy to the device
CopyToDevice(ID, id, N);
//...........................................................................
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// LBM variables
if (rank==0) printf ("Allocating distributions \n");
//......................device distributions.................................
double *f_even,*f_odd;
//...........................................................................
AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size); // Allocate device memory
AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size); // Allocate device memory
//...........................................................................
double *Velocity, *Pressure, *dvcSignDist;
//...........................................................................
AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
//...........................................................................
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages.SDs.data(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
//...........................................................................
// MAIN VARIABLES INITIALIZED HERE
//...........................................................................
//...........................................................................
if (rank==0) printf("Setting the distributions, size = %i\n", N);
//...........................................................................
InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
//......................................................................
//.......................................................................
// 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;
}
}
}
//.......................................................................
if (rank==0 && pBC){
printf("Setting inlet pressure = %f \n", din);
printf("Setting outlet pressure = %f \n", dout);
}
if (pBC && kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
}
if (pBC && kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
}
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,vawx,vawy,vawz,Fo,Re,velocity,err1D,mag_force,vel_prev;
FILE * NONDARCY;
if (rank == 0){
NONDARCY = fopen("nondarcy.csv","a");
fprintf(NONDARCY,"D32 Fx Fy Fz vx vy vz Re Fo\n");
}
Re = 0.f;
// Generate a bunch of points until sufficiently high Re is obtained
while (Re < REYNOLDS_NUMBER){
// Increase the external force and simulate to steady state
Fz = 2.0*Fz;
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 ){
//*************************************************************************
// Fused Color Gradient and Collision
//*************************************************************************
MRT( ID,f_even,f_odd,rlxA,rlxB,Fx,Fy,Fz,Nx,Ny,Nz);
//*************************************************************************
//*************************************************************************
// Pack and send the D3Q19 distributions
ScaLBL_Comm.SendD3Q19(f_even, f_odd);
//*************************************************************************
// Swap the distributions for momentum transport
//*************************************************************************
SwapD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
//*************************************************************************
// Wait for communications to complete and unpack the distributions
ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
//*************************************************************************
if (pBC && kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
}
if (pBC && kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
}
//...................................................................................
DeviceBarrier();
MPI_Barrier(comm);
// Timestep completed!
timestep++;
if (rank==0){
// write out csv file
printf("D32 Fx Fy Fz vx vy vz err1d Fo Re K err\n");
}
if (timestep%500 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
// Way more work than necessary -- this is just to get the solid interfacial area!!
Averages.Initialize();
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
vawx = -Averages.vaw_global(0);
vawy = -Averages.vaw_global(1);
vawz = -Averages.vaw_global(2);
// Compute local measures
err = Re; // previous Reynolds number
D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
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;
err = fabs(Re-err);
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
// *************************************************************************
printf("%f ",D32);
printf("%.5g,%.5g,%.5g ",Fx,Fy,Fz);
printf("%.5g,%.5g,%.5g ",vawx,vawy,vawz);
printf("%.5g ",err1D);
printf("%5g ", Fo);
printf("%.5g ", Re);
printf("%.5g ", Re/Fo);
printf("%.5g\n", err);
}
}
}
// Write steady state variables to csv file
if (rank==0){
fprintf(NONDARCY,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",D32,Fx,Fy,Fz,vawx,vawy,vawz,Re,Fo);
fflush(NONDARCY);
}
}
//************************************************************************/
fclose(NONDARCY);
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(Nx*Ny*Nz)/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);
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();
// ****************************************************
}
}

570
tests/lbpm_segmented_pp.cpp
View File

@ -40,140 +40,146 @@ inline double minmod(double &a, double &b){
inline double Eikonal(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
*
* Reference:
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229
*/
/*
* 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
*
* Reference:
* Min C (2010) On reinitializing level set functions, Journal of Computational Physics 229
*/
int i,j,k;
double dt=0.1;
double Dx,Dy,Dz;
double Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
double sign,norm;
double LocalVar,GlobalVar,LocalMax,GlobalMax;
int i,j,k;
double dt=0.1;
double Dx,Dy,Dz;
double Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
double Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
double sign,norm;
double LocalVar,GlobalVar,LocalMax,GlobalMax;
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 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);
// Arrays to store the second derivatives
DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
// Arrays to store the second derivatives
DoubleArray Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
DoubleArray Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
int count = 0;
while (count < timesteps){
int count = 0;
while (count < timesteps){
// Communicate the halo of values
fillData.fill(Distance);
// Communicate the halo of values
fillData.fill(Distance);
// Compute second order derivatives
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
}
}
}
fillData.fill(Dxx);
fillData.fill(Dyy);
fillData.fill(Dzz);
// Compute second order derivatives
for (k=1;k<Dm.Nz-1;k++){
for (j=1;j<Dm.Ny-1;j++){
for (i=1;i<Dm.Nx-1;i++){
Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
}
}
}
fillData.fill(Dxx);
fillData.fill(Dyy);
fillData.fill(Dzz);
LocalMax=LocalVar=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++){
LocalMax=LocalVar=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;
int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
sign = -1;
if (ID[n] == 1) sign = 1;
sign = -1;
if (ID[n] == 1) sign = 1;
// local second derivative terms
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
// local second derivative terms
Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
/* //............Compute upwind derivatives ...................
/* //............Compute upwind derivatives ...................
Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
*/
Dxp = Distance(i+1,j,k);
Dyp = Distance(i,j+1,k);
Dzp = Distance(i,j,k+1);
*/
Dxp = Distance(i+1,j,k)- Distance(i,j,k) - 0.5*Dxxp;
Dyp = Distance(i,j+1,k)- Distance(i,j,k) - 0.5*Dyyp;
Dzp = Distance(i,j,k+1)- Distance(i,j,k) - 0.5*Dzzp;
Dxm = Distance(i-1,j,k);
Dym = Distance(i,j-1,k);
Dzm = Distance(i,j,k-1);
Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){
if (Dxp > Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
// Compute upwind derivatives for Godunov Hamiltonian
if (sign < 0.0){
if (Dxp + Dxm > 0.f) Dx = Dxp*Dxp;
else Dx = Dxm*Dxm;
if (Dyp > Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dyp + Dym > 0.f) Dy = Dyp*Dyp;
else Dy = Dym*Dym;
if (Dzp > Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
else{
if (Dxp < Dxm) Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
else Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
if (Dzp + Dzm > 0.f) Dz = Dzp*Dzp;
else Dz = Dzm*Dzm;
}
else{
if (Dyp < Dym) Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
else Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
if (Dxp + Dxm < 0.f) Dx = Dxp*Dxp;
else Dx = Dxm*Dxm;
if (Dzp < Dzm) Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
else Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
}
if (Dyp + Dym < 0.f) Dy = Dyp*Dyp;
else Dy = Dym*Dym;
norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
if (norm > 1.0) norm=1.0;
Distance(i,j,k) += dt*sign*(1.0 - norm);
LocalVar += dt*sign*(1.0 - norm);
if (Dzp + Dzm < 0.f) Dz = Dzp*Dzp;
else Dz = Dzm*Dzm;
}
if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
LocalMax = fabs(dt*sign*(1.0 - norm));
}
}
}
//Dx = max(Dxp*Dxp,Dxm*Dxm);
//Dy = max(Dyp*Dyp,Dym*Dym);
//Dz = max(Dzp*Dzp,Dzm*Dzm);
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm);
GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz;
count++;
norm=sqrt(Dx + Dy + Dz);
if (norm > 1.0) norm=1.0;
if (count%50 == 0 && Dm.rank==0 )
printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
Distance(i,j,k) += dt*sign*(1.0 - norm);
LocalVar += dt*sign*(1.0 - norm);
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
LocalMax = fabs(dt*sign*(1.0 - norm));
}
}
}
MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_DOUBLE,MPI_MAX,Dm.Comm);
GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz;
count++;
if (count%50 == 0 && Dm.rank==0 )
printf("Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
if (fabs(GlobalMax) < 1e-5){
if (Dm.rank==0) printf("Exiting with max tolerance of 1e-5 \n");
count=timesteps;
}
}
}
return GlobalVar;
return GlobalVar;
}
@ -182,270 +188,128 @@ int main(int argc, char **argv)
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int BC=0;
char Filename[40];
int xStart,yStart,zStart;
// char fluidValue,solidValue;
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int BC=0;
// char fluidValue,solidValue;
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
ifstream image("Segmented.in");
image >> Filename; // Name of data file containing segmented data
image >> Nx; // size of the binary file
image >> Ny;
image >> Nz;
image >> xStart; // offset for the starting voxel
image >> yStart;
image >> zStart;
}
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);
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
char LocalRankFilename[40];
int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (n=0; n<N; n++) Dm.id[n]=1;
Dm.CommInit(comm);
// Read the phase ID
size_t readID;
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"rb");
readID=fread(Dm.id,1,N,ID);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n");
fclose(ID);
// make sure communication
// Set up layers in x direction
for (k=0; k<nz; k++){
for (j=0; j<ny; j++){
Dm.id[k*nx*ny+j*nx]=1;
Dm.id[k*nx*ny+j*nx+nx-1] = 1;
}
}
for (k=0; k<nz; k++){
for (i=0; i<nx; i++){
Dm.id[k*nx*ny+i]=1;
Dm.id[k*nx*ny+(ny-1)*nx+i] = 1;
}
}
for (j=0; j<ny; j++){
for (i=0; i<nx; i++){
Dm.id[j*nx+i]=1;
Dm.id[nx*ny*(nz-1)+j*nx+i] = 1;
}
}
// Initialize the domain and communication
nx+=2; ny+=2; nz+=2;
int count = 0;
N=nx*ny*nz;
char *id;
id = new char [N];
TwoPhase Averages(Dm);
// DoubleArray Distance(nx,ny,nz);
// DoubleArray Phase(nx,ny,nz);
// Solve for the position of the solid 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;
// Initialize the solid phase
if (Dm.id[n] == 0) id[n] = 0;
else id[n] = 1;
}
}
}
// Initialize the signed distance function
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;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*id[n]-1.0;
}
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);
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
}
}
MeanFilter(Averages.SDs);
double LocalVar, TotalVar;
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx,Dm.Ny*Dm.nprocy),Dm.Nz*Dm.nprocz);
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime);
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
/* // Solve for the position of the non-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;
// Initialize the non-wetting phase
if (Dm.id[n] == 1) id[n] = 1;
else id[n] = 0;
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
}
// Initialize the signed distance function
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;
// Initialize distance to +/- 1
Averages.Phase(i,j,k) = 2.0*id[n]-1.0;
}
}
}
MeanFilter(Averages.Phase);
if (rank==0) printf("Initialized non-wetting phase -- Converting to Signed Distance function \n");
SSO(Averages.Phase,id,Dm,100);
char LocalRankFilename[40];
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;
// Initialize distance to +/- 1
// Dilation of the non-wetting phase
Averages.SDn(i,j,k) = -Averages.Phase(i,j,k);
Averages.Phase(i,j,k) = Averages.SDn(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;
if (Averages.SDs(i,j,k) > 0.0){
if (Averages.Phase(i,j,k) > 0.0){
Dm.id[n] = 2;
}
else{
Dm.id[n] = 1;
}
}
else{
Dm.id[n] = 0;
int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (n=0; n<N; n++) Dm.id[n]=1;
Dm.CommInit(comm);
// Read the phase ID
size_t readID;
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"rb");
readID=fread(Dm.id,1,N,ID);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n");
fclose(ID);
// Initialize the domain and communication
nx+=2; ny+=2; nz+=2;
char *id;
id = new char [N];
TwoPhase Averages(Dm);
// DoubleArray Distance(nx,ny,nz);
// DoubleArray Phase(nx,ny,nz);
int count = 0;
// Solve for the position of the solid 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;
// Initialize the solid phase
if (Dm.id[n] == 0) id[n] = 0;
else id[n] = 1;
}
}
}
// Initialize the signed distance function
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;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*double(id[n])-1.0;
}
}
}
MeanFilter(Averages.SDs);
double LocalVar, TotalVar;
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx,Dm.Ny*Dm.nprocy),Dm.Nz*Dm.nprocz);
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime);
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
if (rank==0) printf("Final variation in signed distance function %f \n",TotalVar);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,N,DIST);
fclose(DIST);
}
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SolidVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "Fluid";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PhaseVar);
SolidVar->name = "Solid";
SolidVar->type = IO::VolumeVariable;
SolidVar->dim = 1;
SolidVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
fillData.copy(Averages.SDn,PhaseVar->data);
fillData.copy(Averages.SDs,SolidVar->data);
fillData.copy(Averages.Label_NWP,BlobIDVar->data);
IO::writeData( 0, meshData, 2, comm );
// sprintf(LocalRankFilename,"Phase.%05i",rank);
// FILE *PHASE = fopen(LocalRankFilename,"wb");
// fwrite(Averages.Phase.get(),8,Averages.Phase.length(),PHASE);
// fclose(PHASE);
double beta = 0.95;
if (rank==0) printf("initializing the system \n");
Averages.UpdateSolid();
Averages.UpdateMeshValues();
Dm.CommunicateMeshHalo(Averages.Phase);
Dm.CommunicateMeshHalo(Averages.SDn);
Dm.CommunicateMeshHalo(Averages.SDs);
int timestep=5;
Averages.Initialize();
if (rank==0) printf("computing phase components \n");
Averages.ComponentAverages();
if (rank==0) printf("sorting phase components \n");
Averages.SortBlobs();
Averages.PrintComponents(timestep);
*/
}
MPI_Barrier(comm);
MPI_Barrier(comm);
MPI_Finalize();
return 0;
return 0;
}

4
workflows/Morphological_Analysis/morphological_analysis_utils.py
View File

@ -1,5 +1,7 @@
#!/usr/bin/env python
# TODO: double check the issue of the "view of the array" for all the following functions
# make sure it is the copy of the array, not the view of the array is used for any
# subprocesses.
import numpy as np
import scipy.stats as stats
import scipy.ndimage.morphology as morphology