diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d561eec4..2e8f8d21 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -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 ) diff --git a/tests/lbpm_morph_pp.cpp b/tests/lbpm_morph_pp.cpp index 27a098a1..2423a53d 100644 --- a/tests/lbpm_morph_pp.cpp +++ b/tests/lbpm_morph_pp.cpp @@ -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 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 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 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 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 +#include +#include +#include +#include +#include +#include + +#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(×tepMax,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 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 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 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(); + // **************************************************** + } +} diff --git a/tests/lbpm_segmented_pp.cpp b/tests/lbpm_segmented_pp.cpp index bc18169e..491a68ea 100644 --- a/tests/lbpm_segmented_pp.cpp +++ b/tests/lbpm_segmented_pp.cpp @@ -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 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 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 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 solidValues; - std::vector nwpValues; - std::string line; + std::vector solidValues; + std::vector 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= 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 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 fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1); - std::vector meshData(1); - meshData[0].meshName = "domain"; - meshData[0].mesh = std::shared_ptr( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) ); - std::shared_ptr PhaseVar( new IO::Variable() ); - std::shared_ptr SolidVar( new IO::Variable() ); - std::shared_ptr 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; } diff --git a/workflows/Morphological_Analysis/morphological_analysis_utils.py b/workflows/Morphological_Analysis/morphological_analysis_utils.py index 9869bb92..8b0eb5e3 100644 --- a/workflows/Morphological_Analysis/morphological_analysis_utils.py +++ b/workflows/Morphological_Analysis/morphological_analysis_utils.py @@ -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