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
2017-07-16 20:13:12 -04:00
parent 63404db74f 37ae5caeac
commit 8dad730f68
14 changed files with 726 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
common/TwoPhase.h
View File

@@ -117,7 +117,7 @@ public:
DoubleArray Gws_global;
//...........................................................................
//...........................................................................
int Nx,Ny,Nz;
int Nx,Ny,Nz;
IntArray PhaseID; // Phase ID array (solid=0, non-wetting=1, wetting=2)
BlobIDArray Label_WP; // Wetting phase label
BlobIDArray Label_NWP; // Non-wetting phase label index (0:nblobs-1)

54
example/Workflow/Eos/eos-Decomp.pbs Normal file
View File

@@ -0,0 +1,54 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N Decomp
#PBS -j oe
#PBS -l walltime=02:30:00,nodes=9
##PBS -l walltime=01:00:00,nodes=18
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=$PROJWORK/geo106/Eos-LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
# directory structure should be set up
# directories must include PATTERN within the name (change as needed)
# each directory must contain
# - Domain.in -- file that specifies parallel domain
# - Segmented.in -- file that specifies input file
# - 8-bit binary file with the digital geometry
LIST=$(ls | grep PATTERN)
# integer labels for the solid phase and the non-wetting phase
SOLID=0
NWP=1
# loop over directories and generate input files for parallel simulation
# performs domain decomposition based on Domain.in
# number of processors must match simulation
NUMPROCS=144
for i in $LIST; do
echo $i
cd $i
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_decomp $SOLID $NWP
mkdir -p MEDIA
cp ID* MEDIA
cd ..
done
exit;

55
example/Workflow/Eos/eos-MorphDrain.pbs Executable file
View File

@@ -0,0 +1,55 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N MorphDrain
#PBS -j oe
##PBS -l walltime=02:00:00,nodes=216
#PBS -l walltime=05:00:00,nodes=18
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/eos/LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
module load python_anaconda
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
cp ../MEDIA/ID* ./
LABEL=$(basename $PWD)
# List of saturation values to generate
SW="0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2"
echo "r sw" > morphdrain.csv
# number of processors to use
NUMPROCS=288
#for r in `cat radius.csv`; do
for sw in $SW; do
echo $sw >> sw.log
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_morphdrain_pp $sw > morphdrain.log
radius=$(grep "Final critical radius" morphdrain.log | sed 's/Final critical radius=//g')
saturation=$(grep "Final saturation" morphdrain.log | sed 's/Final saturation=//g')
echo "$radius $saturation" >> morphdrain.csv
tag="sw_$sw"
DIR=$LABEL"_drain_"$tag
mkdir -p $DIR
cp ID.* $DIR
cp SignDist.* $DIR
cp Domain.in $DIR
cp Color.in $DIR
done
exit;

59
example/Workflow/Eos/eos-MorphOpen.pbs Executable file
View File

@@ -0,0 +1,59 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N MorphOpen
#PBS -j oe
##PBS -l walltime=02:00:00,nodes=216
#PBS -l walltime=08:00:00,nodes=18
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/eos/LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
module load python_anaconda
#module swap cray-mpich2 cray-mpich2/5.6.3
#module load cudatoolkit
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
cp ../MEDIA/ID* ./
LABEL=$(basename $PWD)
rm sw.log
# List of saturations to generate using Morphological opening
SW="0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15"
echo "r sw" > morphopen.csv
NUMPROCS=288
for sw in $SW; do
echo $sw >> sw.log
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_morphopen_pp $sw > morphopen.log
radius=$(grep "Final critical radius" morphopen.log | sed 's/Final critical radius=//g')
saturation=$(grep "Final saturation" morphopen.log | sed 's/Final saturation=//g')
echo "$radius $saturation" >> morphopen.csv
tag="sw_$sw"
DIR=$LABEL"_open_"$tag
mkdir -p $DIR
cp ID.* $DIR
cp SignDist.* $DIR
cp Domain.in $DIR
cp Color.in $DIR
done
exit;

33
example/Workflow/Eos/eos-Preprocess.pbs Normal file
View File

@@ -0,0 +1,33 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N Dist
#PBS -j oe
#PBS -l walltime=02:30:00,nodes=18
##PBS -l walltime=01:00:00,nodes=18
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/eos/LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
# generate distance map from the binary image
# input files are ID.xxxxx
# output files are SignDist.xxxxx
NUMPROCS=288
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_pp
exit;

49
example/Workflow/Eos/eos-Random.pbs Executable file
View File

@@ -0,0 +1,49 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N Random
#PBS -j oe
##PBS -l walltime=02:00:00,nodes=216
#PBS -l walltime=12:00:00,nodes=18
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/eos/LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
module load python_anaconda
#module swap cray-mpich2 cray-mpich2/5.6.3
#module load cudatoolkit
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
LABEL=$(basename $PWD)
NUMPROCS=288
# Generate a bunch of random states
for sat in `seq -w 8 4 92`; do
cp ../MEDIA/ID* ./
sw="0.$sat"
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_random_pp $sw 1
tag="sw_$sw"
DIR=$LABEL"_random_"$tag
mkdir -p $DIR
cp ID.* $DIR
cp SignDist.* $DIR
cp Domain.in $DIR
cp Color.in $DIR
done
exit;

53
example/Workflow/Titan/TitanRelPerm.pbs Normal file
View File

@@ -0,0 +1,53 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N Ryan-Sandstone-RelPerm
#PBS -j oe
#PBS -l walltime=6:00:00,nodes=10660
# Number of nodes should be [# directories] x [GPU per directory]
#
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/titan/LBPM-WIA
source $MODULESHOME/init/bash
module unload PrgEnv-gnu PrgEnv-pgi PrgEnv-cray PrgEnv-intel
module load PrgEnv-gnu
module load cudatoolkit
module load dynamic-link
module load python wraprun
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
export MPICH_RDMA_ENABLED_CUDA=1
export MPICH_MAX_THREAD_SAFETY=multiple
LIST=$(ls | grep "s5_")
NUMSIM=$(ls | grep "s5_" | wc -l)
# Requested number of nodes must be 54 x $NUMSIM !!
echo "Number of simulations to run is $NUMSIM (288 GPU per simulation)"
NUMPROCS=288
# BUIlD THE LAUNCH COMMAND
for dir in $LIST; do
ARGS=$ARGS" : -n $NUMPROCS -N 1 -d 16 --w-cd $dir $LBPM_WIA_INSTALL_DIR/bin/lbpm_color_simulator"
done
echo $ARGS | sed 's/./wraprun/' > launch.log
LAUNCH=$(cat launch.log)
echo "running job "
$LAUNCH
exit;
# aprun -n 512 -N 1 -d 16 $LBPM_WIA_INSTALL_DIR/bin/lbpm_color_simulator >> drainage.log

11
sample_scripts/config_build_eos
View File

@@ -2,7 +2,10 @@
#for var in ${LOADEDMODULES//:/ }; do if [ " ${var///*/}" != " modules" ]; then module unload " ${var///*/}" > /dev/null 2>&1; fi; done
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
module unload PrgEnv-intel
module load PrgEnv-gnu
module unload gcc
module load gcc/5.3.0
module load cray-hdf5-parallel
module load cray-netcdf-hdf5parallel
@@ -18,6 +21,9 @@ module load cmake
export MPICH_RDMA_ENABLED_CUDA=0
echo $GNU_VERSION
module list
# Remove CMake files from previous configures
rm -rf CMake*
@@ -39,9 +45,10 @@ cmake \
-D USE_NETCDF=1 \
-D NETCDF_DIRECTORY=$NETCDF_DIR \
-D HDF5_DIRECTORY=$HDF5_DIR \
-D PREFIX=$MEMBERWORK/geo106/eos-LBPM-WIA \
~/LBPM-WIA
# -D PREFIX=$MEMBERWORK/geo106/eos-LBPM-WIA \
#-D CUDA_HOST_COMPILER="/usr/bin/gcc" \

4
sample_scripts/config_build_titan
View File

@@ -5,6 +5,7 @@ module load cmake
export MPICH_RDMA_ENABLED_CUDA=1
#module swap cray-mpich2 cray-mpich2/5.6.3
module swap PrgEnv-pgi PrgEnv-gnu
module load cray-hdf5 silo
# Remove CMake files from previous configures
rm -rf CMake*
@@ -21,6 +22,9 @@ cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CXX_STD=11 \
-D CUDA_FLAGS="-arch sm_35" \
-D USE_SILO=1 \
-D SILO_DIRECTORY=/sw/xk6/silo/4.8/sles11.1_gnu4.5.3 \
-D HDF5_DIRECTORY=/opt/cray/hdf5/1.8.16/GNU/4.9 \
-D CUDA_HOST_COMPILER="/usr/bin/gcc" \
-D USE_CUDA=1 \
-D USE_TIMER=0 \

2
tests/CMakeLists.txt
View File

@@ -10,8 +10,10 @@ ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphopen_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
ADD_LBPM_EXECUTABLE( lbpm_serial_decomp )
ADD_LBPM_EXECUTABLE( lbpm_disc_pp )
ADD_LBPM_EXECUTABLE( lbpm_captube_pp )
ADD_LBPM_EXECUTABLE( lbpm_plates_pp )
ADD_LBPM_EXECUTABLE( lbpm_squaretube_pp )
ADD_LBPM_EXECUTABLE( lbpm_BlobAnalysis )
ADD_LBPM_EXECUTABLE( TestBubble )

34
tests/lbpm_morphdrain_pp.cpp
View File

@@ -232,10 +232,10 @@ int main(int argc, char **argv)
int Nx = nx;
int Ny = ny;
int Nz = nz;
int GlobalNumber = 1;
double GlobalNumber = 1.f;
int count,countGlobal,totalGlobal;
count = 0;
double count,countGlobal,totalGlobal;
count = 0.f;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
@@ -244,14 +244,14 @@ int main(int argc, char **argv)
else{
// initially saturated with wetting phase
id[n] = 2;
count++;
count+=1.0;
}
}
}
}
// 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));
MPI_Allreduce(&count,&totalGlobal,1,MPI_DOUBLE,MPI_SUM,comm);
double porosity=totalGlobal/(double(nprocx*nprocy*nprocz)*double(nx-2)*double(ny-2)*double(nz-2));
if (rank==0) printf("Media Porosity: %f \n",porosity);
@@ -288,17 +288,17 @@ int main(int argc, char **argv)
while (sw_new > SW && Rcrit_new > 0.99){
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR;// decrease critical radius
Rcrit_new -= deltaR*Rcrit_old;// decrease critical radius
sw_old = sw_new;
sw_diff_old = sw_diff_new;
int Window=round(Rcrit_new);
GlobalNumber = 1;
GlobalNumber = 1.0;
while (GlobalNumber != 0){
while (GlobalNumber > 0){
if (rank==0) printf("GlobalNumber=%i \n",GlobalNumber);
int LocalNumber=GlobalNumber=0;
if (rank==0) printf("GlobalNumber=%f \n",GlobalNumber);
double LocalNumber=GlobalNumber=0.f;
for(k=0; k<Nz; k++){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
@@ -317,7 +317,7 @@ int main(int argc, char **argv)
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_new*Rcrit_new){
LocalNumber++;
LocalNumber+=1.0;
id[nn]=1;
}
}
@@ -407,7 +407,7 @@ int main(int argc, char **argv)
UnpackID(Dm.recvList_YZ, Dm.recvCount_YZ ,recvID_YZ, id);
//......................................................................................
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,comm);
// Layer the inlet with NWP
if (kproc == 0){
@@ -433,19 +433,19 @@ int main(int argc, char **argv)
}
count = 0;
count = 1.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
n=k*Nx*Ny+j*Nx+i;
if (id[n] == 2){
count++;
count+=1.0;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sw_new= double(countGlobal)/totalGlobal;
MPI_Allreduce(&count,&countGlobal,1,MPI_DOUBLE,MPI_SUM,comm);
sw_new= countGlobal/totalGlobal;
sw_diff_new = abs(sw_new-SW);
// if (rank==0){
// printf("Final saturation=%f\n",sw_new);

46
tests/lbpm_morphopen_pp.cpp
View File

@@ -61,6 +61,8 @@ int main(int argc, char **argv)
//printf("Critical radius %f (voxels)\n",Rcrit);
printf("Target saturation %f \n",SW);
}
// }
//.......................................................................
@@ -140,9 +142,9 @@ int main(int argc, char **argv)
if (ReadSignDist != size_t(N)) printf("lbpm_morphdrain_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);
int count,countGlobal,totalGlobal;
count = 0;
double maxdist=0;
double count,countGlobal,totalGlobal;
count = 0.f;
double maxdist=0.f;
double maxdistGlobal;
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
@@ -152,7 +154,7 @@ int main(int argc, char **argv)
else{
// initially saturated with wetting phase
id[n] = 2;
count++;
count+=1.0;
// extract maximum distance for critical radius
if ( SignDist(i,j,k) > maxdist) maxdist=SignDist(i,j,k);
}
@@ -160,9 +162,10 @@ int main(int argc, char **argv)
}
}
// total Global is the number of nodes in the pore-space
MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&count,&totalGlobal,1,MPI_DOUBLE,MPI_SUM,comm);
MPI_Allreduce(&maxdist,&maxdistGlobal,1,MPI_DOUBLE,MPI_MAX,comm);
float porosity=float(totalGlobal)/(nprocx*nprocy*nprocz*(nx-2)*(ny-2)*(nz-2));
double volume=double(nprocx*nprocy*nprocz)*double(nx-2)*double(ny-2)*double(nz-2);
double porosity=totalGlobal/volume;
if (rank==0) printf("Media Porosity: %f \n",porosity);
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);\
@@ -342,20 +345,24 @@ int main(int argc, char **argv)
double deltaR=0.05; // amount to change the radius in voxel units
double Rcrit_old;
int GlobalNumber = 1;
double GlobalNumber = 1.f;
int imin,jmin,kmin,imax,jmax,kmax;
Rcrit_new = maxdistGlobal;
if (argc>2){
Rcrit_new = strtod(argv[2],NULL);
if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);
}
while (sw_new > SW)
{
sw_diff_old = sw_diff_new;
sw_old = sw_new;
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR;
Rcrit_new -= deltaR*Rcrit_old;
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;
double LocalNumber=0.f;
for(k=0; k<Nz; k++){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){
@@ -374,7 +381,7 @@ int main(int argc, char **argv)
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_new*Rcrit_new){
LocalNumber++;
LocalNumber+=1.0;
id[nn]=1;
}
}
@@ -465,28 +472,27 @@ int main(int argc, char **argv)
UnpackID(Dm.recvList_YZ, Dm.recvCount_YZ ,recvID_YZ, id);
//......................................................................................
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_INT,MPI_SUM,comm);
MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,comm);
count = 0;
count = 0.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
n=k*Nx*Ny+j*Nx+i;
if (id[n] == 2){
count++;
count+=1.0;
}
}
}
}
MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,comm);
sw_new = float(countGlobal)/totalGlobal;
MPI_Allreduce(&count,&countGlobal,1,MPI_DOUBLE,MPI_SUM,comm);
sw_new = 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 (rank==0){
printf("Final saturation=%f\n",sw_new);
printf("Final critical radius=%f\n",Rcrit_new);
}
}
if (sw_diff_new<sw_diff_old){

204
tests/lbpm_plates_pp.cpp Normal file
View File

@@ -0,0 +1,204 @@
#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"
int main(int argc, char **argv)
{
//*****************************************
// ***** 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;
int sendtag,recvtag;
//*****************************************
// 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 TubeRadius =15.0;
double Width;
TubeRadius=strtod(argv[1],NULL);
WIDTH=strtod(argv[2],NULL);
int BC = 0;
if (rank == 0){
printf("********************************************************\n");
printf("Generate 3D parallel plates geometry with radius = %f voxels \n",TubeRadius);
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
int i,j,k,n;
// pmmc threshold values
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);
printf("nprocz = %i \n",nprocz);
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
if (rank==0){
printf("********************************************************\n");
printf("Sub-domain size = %i x %i x %i\n",Nz,Nz,Nz);
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
printf("********************************************************\n");
}
// Initialized domain and averaging framework for Two-Phase Flow
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
Dm.CommInit(comm);
TwoPhase Averages(Dm);
MPI_Barrier(comm);
Nz += 2;
Nx = Ny = Nz; // Cubic domain
int N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
//.......................................................................
// 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 pore_vol;
sum=0;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// Cylindrical capillary tube aligned with the z direction
Averages.SDs(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2));
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Averages.SDs(i,j,k) < WIDTH){
id[n] = 2;
sum++;
}
else {
id[n] = 1;
sum++;
}
}
}
}
// Compute the pore volume
sum_local = 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 (id[n] > 0){
sum_local += 1.0;
}
}
}
}
MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
//.........................................................
// 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;
//.........................................................
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);
fclose(ID);
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();
// ****************************************************
}

160
tests/lbpm_serial_decomp.cpp Normal file
View File

@@ -0,0 +1,160 @@
/*
* Pre-processor to generate signed distance function from segmented data
* segmented data should be stored in a raw binary file as 1-byte integer (type char)
* will output distance functions for phases
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include "common/Array.h"
#include "common/Domain.h"
int main(int argc, char **argv)
{
int rank=0;
bool MULTINPUT=false;
int NWP,SOLID,rank_offset;
SOLID=atoi(argv[1]);
NWP=atoi(argv[2]);
if (rank==0){
printf("Solid Label: %i \n",SOLID);
printf("NWP Label: %i \n",NWP);
}
if (argc > 3){
rank_offset = atoi(argv[3]);
}
else{
MULTINPUT=true;
rank_offset=0;
}
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
uint64_t Nx,Ny,Nz;
uint64_t i,j,k,n;
int BC=0;
char Filename[40];
uint64_t xStart,yStart,zStart;
// char fluidValue,solidValue;
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;
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;
}
nprocs=nprocx*nprocy*nprocz;
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
uint64_t SIZE = Nx*Ny*Nz;
SegData = new char[SIZE];
FILE *SEGDAT = fopen(Filename,"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename);
}
// Get the rank info
uint64_t N = (nx+2)*(ny+2)*(nz+2);
// number of sites to use for periodic boundary condition transition zone
uint64_t z_transition_size = (nprocz*nz - (Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0;
char LocalRankFilename[40];
char *loc_id;
loc_id = new char [(nx+2)*(ny+2)*(nz+2)];
// Set up the sub-domains
if (rank==0){
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
printf("Size of transition region: %i \n", z_transition_size);
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*nprocx*nprocy + jp*nprocx + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
uint64_t x = xStart + ip*nx + i-1;
uint64_t y = yStart + jp*ny + j-1;
// uint64_t z = zStart + kp*nz + k-1;
uint64_t z = zStart + kp*nz + k-1 - z_transition_size;
if (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
uint64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
uint64_t nglobal = z*Nx*Ny+y*Nx+x;
loc_id[nlocal] = SegData[nglobal];
}
}
}
// relabel the data
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
n = k*(nx+2)*(ny+2) + j*(nx+2) + i;;
if (loc_id[n]==char(SOLID)) loc_id[n] = 0;
else if (loc_id[n]==char(NWP)) loc_id[n] = 1;
else loc_id[n] = 2;
}
}
}
// Write the data for this rank data
sprintf(LocalRankFilename,"ID.%05i",rnk+rank_offset);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(loc_id,1,(nx+2)*(ny+2)*(nz+2),ID);
fclose(ID);
}
}
}
}
printf("Domain decomposition completed successfully \n");
return 0;
}