Adding a nonnewtonian example and test
This commit is contained in:
cpf
@@ -1,4 +1,5 @@
|
||||
# Copy the examples to the install folder
|
||||
INSTALL_EXAMPLE (NonNewtonianChannelFlow )
|
||||
INSTALL_EXAMPLE( Bubble )
|
||||
INSTALL_EXAMPLE( ConstrainedBubble )
|
||||
INSTALL_EXAMPLE( Piston )
|
||||
|
||||
@@ -0,0 +1,4 @@
|
||||
1 1 1
|
||||
100 100 100
|
||||
229
|
||||
1.0 1.0 1.0
|
||||
@@ -0,0 +1,4 @@
|
||||
1.0
|
||||
0.0 0.0 1.0e-5
|
||||
0 0 1.0 1.0
|
||||
1000 1000 1.0e-5
|
||||
@@ -0,0 +1,11 @@
|
||||
# 1. Edit the Domain.in file based on the desired system size
|
||||
|
||||
|
||||
# 2. Run the pre-processor
|
||||
export TUBEWIDTH=20
|
||||
export LAYERWIDTH=0
|
||||
mpirun -np 1 ../../tests/lbpm_plates_pp $TUBEWIDTH $LAYERWIDTH
|
||||
|
||||
|
||||
# 3. Run single-phase simulation within the domain
|
||||
mpirun -np 1 ../../tests/lbpm_permeability_simulator
|
||||
@@ -1,4 +1,4 @@
|
||||
# Configure
|
||||
#Configure
|
||||
|
||||
rm -rf Cmake*
|
||||
|
||||
@@ -7,16 +7,20 @@ cmake \
|
||||
-D CMAKE_ECLIPSE_GENERATE_SOURCE_PROJECT=TRUE \
|
||||
-D CMAKE_ECLIPSE_MAKE_ARGUMENTS=-j8 \
|
||||
-D CMAKE_C_COMPILER:PATH=/home/christopher/openmpi/install_dir/bin/mpicc \
|
||||
-D CMAKE_CXX_COMPILER:PATH=/home/christopher/openmpi/install_dir/bin/mpicxx \
|
||||
i -D CMAKE_CXX_COMPILER:PATH=/home/christopher/openmpi/install_dir/bin/mpicxx \
|
||||
-D USE_MPI=1 \
|
||||
-D USE_DOXYGEN=false \
|
||||
-D CMAKE_C_FLAGS="-std=gnu++11 -w" \
|
||||
-D CMAKE_CXX_FLAGS="-std=gnu++11 -w" \
|
||||
# -D CMAKE_C_FLAGS="-std=gnu++11 -w" \
|
||||
# -D CMAKE_CXX_FLAGS="-std=gnu++11 -w" \
|
||||
-D MPI_COMPILER:BOOL=TRUE \
|
||||
-D MPIEXEC=/home/christopher/openmpi/install_dir/bin/mpirun \
|
||||
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
|
||||
-D CMAKE_BUILD_TYPE:STRING=Release \
|
||||
-D CXX_STD=11 \
|
||||
-D HDF5_DIRECTORY="/usr/lib/x86_64-linux-gnu/hdf5/openmpi/" \
|
||||
-D USE_SILO=1 \
|
||||
-D SILO_DIRECTORY="/home/christopher/silo-4.10.2/" \
|
||||
-D SILO_LIB="/home/christopher/silo-4.10.2/lib/libsiloh5.a" \
|
||||
-D CUDA_FLAGS="-arch sm_60 -Xcompiler -fPIC -std=c++11" \
|
||||
-D CUDA_HOST_COMPILER="/usr/bin/gcc-4.9" \
|
||||
-D USE_CUDA=1 \
|
||||
@@ -26,4 +30,4 @@ cmake \
|
||||
/home/christopher/repos/LBPM-WIA
|
||||
|
||||
make -j8 install
|
||||
make visit
|
||||
#make visit
|
||||
|
||||
@@ -1,4 +1,5 @@
|
||||
# Copy files for the tests
|
||||
ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator )
|
||||
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
|
||||
ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
|
||||
ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
|
||||
|
||||
@@ -0,0 +1,977 @@
|
||||
#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"
|
||||
#include "ProfilerApp.h"
|
||||
#include "threadpool/thread_pool.h"
|
||||
|
||||
#include "lbpm_nonnewtonian_simulator.h"
|
||||
|
||||
//#define WRITE_SURFACES
|
||||
|
||||
/*
|
||||
* Simulator for single-phase non-newtonian flow
|
||||
* James E. McClure 2013-2014 & Christopher P. Fowler 2017
|
||||
*/
|
||||
|
||||
using namespace std;
|
||||
|
||||
//*************************************************************************
|
||||
// Implementation of Steady State Single-Phase LBM for permeability measurement
|
||||
//*************************************************************************
|
||||
inline void PackID(int *list, int count, char *sendbuf, char *ID){
|
||||
// Fill in the phase ID values from neighboring processors
|
||||
// This packs up the values that need to be sent from one processor to another
|
||||
int idx,n;
|
||||
|
||||
for (idx=0; idx<count; idx++){
|
||||
n = list[idx];
|
||||
sendbuf[idx] = ID[n];
|
||||
}
|
||||
}
|
||||
//***************************************************************************************
|
||||
|
||||
inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
|
||||
// Fill in the phase ID values from neighboring processors
|
||||
// This unpacks the values once they have been recieved from neighbors
|
||||
int idx,n;
|
||||
|
||||
for (idx=0; idx<count; idx++){
|
||||
n = list[idx];
|
||||
ID[n] = recvbuf[idx];
|
||||
}
|
||||
}
|
||||
|
||||
//***************************************************************************************
|
||||
|
||||
inline void ZeroHalo(double *Data, int Nx, int Ny, int Nz)
|
||||
{
|
||||
int i,j,k,n;
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
i=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
i=Nx-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
for (k=0;k<Nz;k++){
|
||||
for (i=0;i<Nx;i++){
|
||||
j=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
j=Ny-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
k=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
k=Nz-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
//***************************************************************************************
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
//*****************************************
|
||||
// ***** MPI STUFF ****************
|
||||
//*****************************************
|
||||
// Initialize MPI
|
||||
//MPI_Init(&argc,&argv);
|
||||
|
||||
/*
|
||||
* Definitely seems to be an issue - let's hope James gets back to me...
|
||||
*/
|
||||
int provided_thread_support = -1;
|
||||
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
|
||||
MPI_Comm comm;
|
||||
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
|
||||
int rank = comm_rank(comm);
|
||||
int nprocs = comm_size(comm);
|
||||
|
||||
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
|
||||
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
|
||||
|
||||
|
||||
|
||||
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
|
||||
|
||||
// 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];
|
||||
|
||||
if (rank == 0){
|
||||
printf("********************************************************\n");
|
||||
printf("Running Single Phase Non-Newtonian Calculation \n");
|
||||
printf("********************************************************\n");
|
||||
}
|
||||
|
||||
// Variables that specify the computational domain
|
||||
string FILENAME;
|
||||
int Nx,Ny,Nz; // local sub-domain size
|
||||
int nspheres; // number of spheres in the packing
|
||||
double Lx,Ly,Lz; // Domain length
|
||||
double D = 1.0; // reference length for non-dimensionalization
|
||||
// Color Model parameters
|
||||
int timestepMax, interval;
|
||||
double tau,Fx,Fy,Fz,tol,err;
|
||||
double din,dout;
|
||||
bool pBC,Restart;
|
||||
int i,j,k,n;
|
||||
|
||||
|
||||
/*
|
||||
* Analysis flags
|
||||
*/
|
||||
int RESTART_INTERVAL=20000;
|
||||
int BLOB_ANALYSIS_INTERVAL=1000;
|
||||
int timestep = -1;
|
||||
|
||||
/*
|
||||
* Read file data in
|
||||
*
|
||||
*/
|
||||
|
||||
if (rank==0){
|
||||
//.............................................................
|
||||
// READ SIMULATION PARMAETERS FROM INPUT FILE
|
||||
//.............................................................
|
||||
ifstream input("Permeability.in");
|
||||
// Line 1: model parameters (tau, alpha, beta, das, dbs)
|
||||
input >> tau; // Viscosity parameter
|
||||
// Line 2: External force components (Fx,Fy, Fz)
|
||||
input >> Fx;
|
||||
input >> Fy;
|
||||
input >> Fz;
|
||||
// Line 3: Pressure Boundary conditions
|
||||
input >> Restart;
|
||||
input >> pBC;
|
||||
input >> din;
|
||||
input >> dout;
|
||||
// Line 4: time-stepping criteria
|
||||
input >> timestepMax; // max no. of timesteps
|
||||
input >> interval; // restart interval
|
||||
input >> tol; // error tolerance
|
||||
//.............................................................
|
||||
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> Nx;
|
||||
domain >> Ny;
|
||||
domain >> Nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
//.......................................................................
|
||||
|
||||
}
|
||||
|
||||
/*
|
||||
* Broadcast
|
||||
*/
|
||||
|
||||
// **************************************************************
|
||||
// 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;
|
||||
// **************************************************************
|
||||
// **************************************************************
|
||||
|
||||
/*
|
||||
* Set up rank info struct
|
||||
*/
|
||||
|
||||
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
|
||||
|
||||
MPI_Barrier(comm);
|
||||
|
||||
/*
|
||||
* Set up the relaxation rates and STATIC VISCOSITY
|
||||
*/
|
||||
|
||||
double rlxA = 1.f/tau;
|
||||
double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
|
||||
double viscosity=(tau-0.5)/3.0;
|
||||
|
||||
/*
|
||||
* Debug block 1
|
||||
*/
|
||||
printf("\npBC=%d (an int) \n",pBC);
|
||||
printf("viscosity=%f\n",viscosity);
|
||||
|
||||
/*
|
||||
* Check processor counts
|
||||
*/
|
||||
|
||||
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!");
|
||||
}
|
||||
|
||||
/*
|
||||
* Display what we've got thus far
|
||||
*/
|
||||
|
||||
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");
|
||||
}
|
||||
|
||||
/*
|
||||
* Initialized domain and averaging framework for Two-Phase flow
|
||||
*/
|
||||
|
||||
// not needed right now
|
||||
|
||||
// Initialized domain and averaging framework for Two-Phase Flow
|
||||
int BC=pBC;
|
||||
|
||||
printf("BC=pBC=%d (an int)\n",BC);
|
||||
|
||||
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC); /* 1 */
|
||||
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
|
||||
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
|
||||
Dm.CommInit(comm); /* 2 */
|
||||
|
||||
Domain Mask(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);
|
||||
|
||||
/*
|
||||
* Debug block 2
|
||||
*/
|
||||
printf("LocalRankString=%s\n",LocalRankString);
|
||||
printf("LocalRestartFile=%s\n",LocalRestartFile);
|
||||
printf("LocalRankFilename=%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) {
|
||||
printf("tripped if (pBC) statement\n");
|
||||
iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
|
||||
}
|
||||
double porosity, pore_vol;
|
||||
|
||||
|
||||
/*
|
||||
* Debug block 3
|
||||
*/
|
||||
printf("iVol_global=%f\n",iVol_global);
|
||||
|
||||
//...........................................................................
|
||||
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; /* 3 */
|
||||
|
||||
//.......................................................................
|
||||
// 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++;
|
||||
}
|
||||
}
|
||||
}
|
||||
} /* 4 */
|
||||
|
||||
/*
|
||||
* Initialize from segmented data 5
|
||||
*
|
||||
*/
|
||||
|
||||
// not needed at the moment
|
||||
|
||||
|
||||
/*
|
||||
* Debug block 4
|
||||
*/
|
||||
printf("sum=%d\n",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); /* 6 */
|
||||
//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 external boundary conditions are applied remove solid
|
||||
if (pBC && kproc == 0){
|
||||
printf("Tripped if (pcB && kproc == 0)\n");
|
||||
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){
|
||||
printf("Tripped if (pcB && kproc == nprocz-1)\n");
|
||||
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; /* 7 */
|
||||
//.........................................................
|
||||
|
||||
/*
|
||||
* To use a mask or not - that is the question!
|
||||
*/
|
||||
// maybe mask? /* 8 */
|
||||
|
||||
// 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);
|
||||
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
|
||||
Mask.CommInit(comm);
|
||||
|
||||
//...........................................................................
|
||||
if (rank==0) printf ("Create ScaLBL_Communicator \n");
|
||||
// Create a communicator for the device
|
||||
// ScaLBL_Communicator ScaLBL_Comm(Dm);
|
||||
ScaLBL_Communicator ScaLBL_Comm(Mask); /* 9 */
|
||||
|
||||
// set reservoirs (not needed, right?)
|
||||
if (pBC > 0){
|
||||
for ( k=0;k<Nz;k++){
|
||||
for ( j=0;j<Ny;j++){
|
||||
for ( i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
if (Dm.kproc==0 && k==0) id[n]=1;
|
||||
if (Dm.kproc==0 && k==1) id[n]=1;
|
||||
if (Dm.kproc==nprocz-1 && k==Nz-2) id[n]=2;
|
||||
if (Dm.kproc==nprocz-1 && k==Nz-1) id[n]=2;
|
||||
Mask.id[n] = id[n];
|
||||
}
|
||||
}
|
||||
}
|
||||
} /* 10 */
|
||||
|
||||
//...........device phase ID.................................................
|
||||
if (rank==0) printf ("Copying phase ID to device \n");
|
||||
char *ID;
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
if (i==0 || i==Nx-1 || j==0 || j==Ny-1 || k==0 || k==Nz-1) id[n] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Copy to the device
|
||||
ScaLBL_CopyToDevice(ID, id, N);
|
||||
ScaLBL_DeviceBarrier();
|
||||
//........................................................................... 11
|
||||
//...........................................................................
|
||||
// MAIN VARIABLES ALLOCATED HERE
|
||||
//...........................................................................
|
||||
// LBM variables
|
||||
if (rank==0) printf ("Allocating distributions \n");
|
||||
//......................device distributions.................................
|
||||
double *f_even,*f_odd;
|
||||
//...........................................................................
|
||||
ScaLBL_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size); // Allocate device memory
|
||||
//...........................................................................
|
||||
double *Velocity, *Pressure, *dvcSignDist;
|
||||
//...........................................................................
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
|
||||
//...........................................................................
|
||||
|
||||
// Copy signed distance for device initialization
|
||||
ScaLBL_CopyToDevice(dvcSignDist, Averages->SDs.data(), dist_mem_size);
|
||||
//...........................................................................
|
||||
|
||||
int logcount = 0; // number of surface write-outs /* 12 */
|
||||
|
||||
/*
|
||||
* Display simulation metrics:
|
||||
*/
|
||||
if (rank == 0) {
|
||||
printf("Displaying simulation metrics... \n");
|
||||
}
|
||||
|
||||
|
||||
//...........................................................................
|
||||
// MAIN VARIABLES INITIALIZED HERE
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
if (rank==0) printf("Setting the distributions, size = %i\n", N);
|
||||
//...........................................................................
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q19_Init(ID, f_even, f_odd, Nx, Ny, Nz);
|
||||
ScaLBL_DeviceBarrier();
|
||||
//...................................................................... /* 13 */
|
||||
|
||||
|
||||
if (Restart == true){
|
||||
|
||||
|
||||
if (rank==0){
|
||||
printf("Reading restart file! \n");
|
||||
ifstream restart("Restart.txt");
|
||||
if (restart.is_open()){
|
||||
restart >> timestep;
|
||||
printf("Restarting from timestep =%i \n",timestep);
|
||||
}
|
||||
else{
|
||||
printf("WARNING:No Restart.txt file, setting timestep=0 \n");
|
||||
timestep=5;
|
||||
}
|
||||
}
|
||||
MPI_Bcast(×tep,1,MPI_INT,0,comm);
|
||||
|
||||
// Read in the restart file to CPU buffers
|
||||
double *cDen = new double[2*N];
|
||||
double *cDistEven = new double[10*N];
|
||||
double *cDistOdd = new double[9*N];
|
||||
ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
|
||||
// Copy the restart data to the GPU
|
||||
ScaLBL_CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
|
||||
ScaLBL_CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
|
||||
// ScaLBL_CopyToDevice(Den,cDen,2*N*sizeof(double)); /* Two-phase stuff */
|
||||
ScaLBL_DeviceBarrier();
|
||||
delete [] cDen;
|
||||
delete [] cDistEven;
|
||||
delete [] cDistOdd;
|
||||
MPI_Barrier(comm);
|
||||
} /* 14 */
|
||||
|
||||
// //......................................................................
|
||||
// ScaLBL_D3Q7_Init(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
|
||||
// ScaLBL_D3Q7_Init(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
|
||||
// ScaLBL_DeviceBarrier();
|
||||
// MPI_Barrier(comm); /* 15 */
|
||||
|
||||
//.......................................................................
|
||||
// Once phase has been initialized, map solid to account for 'smeared' interface
|
||||
//for (i=0; i<N; i++) Averages->SDs(i) -= (1.0);
|
||||
// Make sure the id match for the two domains
|
||||
for (i=0; i<N; i++) Dm.id[i] = Mask.id[i];
|
||||
//....................................................................... /* 16 */
|
||||
|
||||
//.......................................................................
|
||||
// Finalize setup for averaging domain
|
||||
Averages->UpdateSolid(); /* 17 */
|
||||
|
||||
|
||||
|
||||
// //.......................................................................
|
||||
//
|
||||
// //*************************************************************************
|
||||
// // Compute the phase indicator field and reset Copy, Den
|
||||
// //*************************************************************************
|
||||
// ScaLBL_ComputePhaseField(ID, Phi, Den, N);
|
||||
// //*************************************************************************
|
||||
// ScaLBL_DeviceBarrier();
|
||||
// ScaLBL_Comm.SendHalo(Phi);
|
||||
// ScaLBL_Comm.RecvHalo(Phi);
|
||||
// ScaLBL_DeviceBarrier();
|
||||
// MPI_Barrier(comm);
|
||||
// //************************************************************************* /* 18 */
|
||||
|
||||
|
||||
if (rank==0 && pBC){
|
||||
printf("Setting inlet pressure = %f \n", din);
|
||||
printf("Setting outlet pressure = %f \n", dout);
|
||||
}
|
||||
if (pBC && kproc == 0) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
}
|
||||
|
||||
if (pBC && kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
} /* 19 */
|
||||
|
||||
// will have to fill in stuff for 2 phase later
|
||||
|
||||
//timestepMax = 500;
|
||||
|
||||
if (rank==0) printf("********************************************************\n");
|
||||
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
|
||||
|
||||
//...........................................................................
|
||||
// Copy the data for for the analysis timestep
|
||||
//...........................................................................
|
||||
// Copy the phase from the GPU -> CPU
|
||||
//...........................................................................
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
// ScaLBL_CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Press.data(),Pressure,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_x.data(),&Velocity[0],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_y.data(),&Velocity[N],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],N*sizeof(double));
|
||||
//........................................................................... /* 20 */
|
||||
|
||||
//.......create and start timer............
|
||||
double starttime,stoptime,cputime;
|
||||
MPI_Barrier(comm);
|
||||
starttime = MPI_Wtime();
|
||||
|
||||
/*
|
||||
* Create the thread pool
|
||||
*
|
||||
*/
|
||||
|
||||
// Create the thread pool
|
||||
int N_threads = 4;
|
||||
if ( provided_thread_support < MPI_THREAD_MULTIPLE )
|
||||
N_threads = 0;
|
||||
if ( N_threads > 0 ) {
|
||||
// Set the affinity
|
||||
int N_procs = ThreadPool::getNumberOfProcessors();
|
||||
std::vector<int> procs(N_procs);
|
||||
for (int i=0; i<N_procs; i++)
|
||||
procs[i] = i;
|
||||
ThreadPool::setProcessAffinity(procs);
|
||||
}
|
||||
ThreadPool tpool(N_threads); /* 21 */
|
||||
|
||||
printf("N_threads=%d\n",N_threads);
|
||||
|
||||
/*
|
||||
* Create the MeshDataStruct - seems particularly important
|
||||
*/
|
||||
|
||||
// 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> PressVar( new IO::Variable() );
|
||||
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
|
||||
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
|
||||
PhaseVar->name = "phase";
|
||||
PhaseVar->type = IO::VariableType::VolumeVariable;
|
||||
PhaseVar->dim = 1;
|
||||
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
|
||||
meshData[0].vars.push_back(PhaseVar);
|
||||
PressVar->name = "Pressure";
|
||||
PressVar->type = IO::VariableType::VolumeVariable;
|
||||
PressVar->dim = 1;
|
||||
PressVar->data.resize(Nx-2,Ny-2,Nz-2);
|
||||
meshData[0].vars.push_back(PressVar);
|
||||
SignDistVar->name = "SignDist";
|
||||
SignDistVar->type = IO::VariableType::VolumeVariable;
|
||||
SignDistVar->dim = 1;
|
||||
SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
|
||||
meshData[0].vars.push_back(SignDistVar);
|
||||
BlobIDVar->name = "BlobID";
|
||||
BlobIDVar->type = IO::VariableType::VolumeVariable;
|
||||
BlobIDVar->dim = 1;
|
||||
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
|
||||
meshData[0].vars.push_back(BlobIDVar); /* 22 */
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//.........................................
|
||||
|
||||
double D32,Fo,Re,velocity,err1D,mag_force,vel_prev;
|
||||
err = vel_prev = 1.0;
|
||||
if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
|
||||
//************ MAIN ITERATION LOOP ***************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
BlobIDstruct last_ids, last_index;
|
||||
BlobIDList last_id_map;
|
||||
writeIDMap(ID_map_struct(),0,id_map_filename);
|
||||
AnalysisWaitIdStruct work_ids;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
double beta = 0;
|
||||
// double * Phi;
|
||||
// double * Den;
|
||||
|
||||
|
||||
|
||||
|
||||
while (timestep < timestepMax && err > tol ){
|
||||
|
||||
//*************************************************************************
|
||||
// Fused Color Gradient and Collision
|
||||
//*************************************************************************
|
||||
ScaLBL_D3Q19_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
|
||||
//*************************************************************************
|
||||
ScaLBL_D3Q19_Swap(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) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
}
|
||||
|
||||
if (pBC && kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
}
|
||||
//...................................................................................
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
|
||||
// Timestep completed!
|
||||
timestep++;
|
||||
|
||||
/// Perform the analysis
|
||||
run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
|
||||
Nx,Ny,Nz,pBC,err,Pressure,Velocity,ID,f_even,f_odd,
|
||||
LocalRestartFile,meshData,fillData,tpool,work_ids);
|
||||
|
||||
|
||||
}
|
||||
//************************************************************************/
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
stoptime = MPI_Wtime();
|
||||
if (rank==0) printf("-------------------------------------------------------------------\n");
|
||||
// Compute the walltime per timestep
|
||||
cputime = (stoptime - starttime)/timestep;
|
||||
// Performance obtained from each node
|
||||
double MLUPS = double(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();
|
||||
//****************************************************
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// Scrap
|
||||
|
||||
// if (rank==0){
|
||||
// ************* DIMENSIONLESS FORCHEIMER EQUATION *************************
|
||||
// Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller
|
||||
// Description of Non-Darcy Flows in Porous Medium Systems
|
||||
// Physical Review E 87 (3), 033012
|
||||
// Fo := density*D32^3*(density*force) / (viscosity^2)
|
||||
// Re := density*D32*velocity / viscosity
|
||||
// Fo = a*Re + b*Re^2
|
||||
// *************************************************************************
|
||||
//viscosity = (tau-0.5)*0.333333333333333333;
|
||||
/*
|
||||
* Original formula for D32: D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
|
||||
*/
|
||||
// D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
|
||||
// printf("Dm.Volume=%f Averages.vol_w_global=%f Averages.As_global=%f \n",Dm.Volume,Averages.vol_w_global,Averages.As_global);
|
||||
// D32 = 6.0*(Dm.Volume-Averages.vol_w_global);
|
||||
// printf("Sauter Mean Diameter = %f \n",D32);
|
||||
// mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
|
||||
// Fo = D32*D32*D32*mag_force/viscosity/viscosity;
|
||||
// .... 1-D flow should be aligned with force ...
|
||||
// velocity = vawx*Fx/mag_force + vawy*Fy/mag_force + vawz*Fz/mag_force;
|
||||
// err1D = fabs(velocity-sqrt(vawx*vawx+vawy*vawy+vawz*vawz))/velocity;
|
||||
//.......... Computation of the Reynolds number Re ..............
|
||||
// Re = D32*velocity/viscosity;
|
||||
// printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
|
||||
// printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
|
||||
// printf("Relative error for 1D representation: %.5g \n",err1D);
|
||||
// printf("Dimensionless force: %5g \n", Fo);
|
||||
// printf("Reynolds number: %.5g \n", Re);
|
||||
// printf("Dimensionless Permeability (k/D^2): %.5g \n", Re/Fo);
|
||||
// }
|
||||
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* if (timestep%5 == 0){
|
||||
//...........................................................................
|
||||
// Copy the data for for the analysis timestep
|
||||
//...........................................................................
|
||||
// Copy the phase from the GPU -> CPU
|
||||
//...........................................................................
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
|
||||
ScaLBL_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();
|
||||
|
||||
double vawx = -Averages.vaw_global(0);
|
||||
double vawy = -Averages.vaw_global(1);
|
||||
double vawz = -Averages.vaw_global(2);
|
||||
|
||||
|
||||
if (rank==0){
|
||||
mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
|
||||
// .... 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;
|
||||
//printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
|
||||
printf("vel_z=%.5g\n",vawz);
|
||||
//printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
|
||||
printf("Relative error for 1D representation: %.5g \n",err1D);
|
||||
}
|
||||
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
|
||||
// // 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;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// //.......................................................................
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
//
|
||||
@@ -0,0 +1,486 @@
|
||||
// Run the analysis, blob identification, and write restart files
|
||||
#include "common/Array.h"
|
||||
#include "common/Communication.h"
|
||||
#include "common/MPI_Helpers.h"
|
||||
#include "IO/MeshDatabase.h"
|
||||
|
||||
//#define ANALYSIS_INTERVAL 6
|
||||
#define ANALYSIS_INTERVAL 1000
|
||||
#define BLOBID_INTERVAL 1000
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
|
||||
CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20 };
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
void DeleteArray( const TYPE *p )
|
||||
{
|
||||
delete [] p;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
// Structure used to store ids
|
||||
struct AnalysisWaitIdStruct {
|
||||
ThreadPool::thread_id_t blobID;
|
||||
ThreadPool::thread_id_t analysis;
|
||||
ThreadPool::thread_id_t vis;
|
||||
ThreadPool::thread_id_t restart;
|
||||
};
|
||||
|
||||
|
||||
|
||||
// Helper class to write the restart file from a seperate thread
|
||||
class WriteRestartWorkItem: public ThreadPool::WorkItem
|
||||
{
|
||||
public:
|
||||
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_,
|
||||
std::shared_ptr<double> cDistEven_, std::shared_ptr<double>cDistOdd_, int N_ ):
|
||||
filename(filename_), cDen(cDen_), cDistEven(cDistEven_), cDistOdd(cDistOdd_), N(N_) {}
|
||||
virtual void run() {
|
||||
PROFILE_START("Save Checkpoint",1);
|
||||
WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
|
||||
PROFILE_STOP("Save Checkpoint",1);
|
||||
};
|
||||
virtual bool has_result() const { return false; }
|
||||
private:
|
||||
WriteRestartWorkItem();
|
||||
const char* filename;
|
||||
std::shared_ptr<double> cDen, cDistEven, cDistOdd;
|
||||
const int N;
|
||||
};
|
||||
|
||||
|
||||
// Helper class to compute the blob ids
|
||||
static const std::string id_map_filename = "lbpm_id_map.txt";
|
||||
typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
|
||||
typedef std::shared_ptr<std::vector<BlobIDType> > BlobIDList;
|
||||
|
||||
|
||||
|
||||
//
|
||||
//class BlobIdentificationWorkItem1: public ThreadPool::WorkItem
|
||||
//{
|
||||
//public:
|
||||
// BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
||||
// std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
||||
// BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
|
||||
// timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
||||
// phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
|
||||
// {
|
||||
// MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
|
||||
// }
|
||||
// ~BlobIdentificationWorkItem1() { MPI_Comm_free(&newcomm); }
|
||||
// virtual void run() {
|
||||
// // Compute the global blob id and compare to the previous version
|
||||
// PROFILE_START("Identify blobs",1);
|
||||
// double vF = 0.0;
|
||||
// double vS = -1.0; // one voxel buffer region around solid
|
||||
// IntArray& ids = new_index->second;
|
||||
// new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,newcomm);
|
||||
// PROFILE_STOP("Identify blobs",1);
|
||||
// }
|
||||
// virtual bool has_result() const { return false; }
|
||||
//private:
|
||||
// BlobIdentificationWorkItem1();
|
||||
// int timestep;
|
||||
// int Nx, Ny, Nz;
|
||||
// const RankInfoStruct& rank_info;
|
||||
// std::shared_ptr<const DoubleArray> phase;
|
||||
// const DoubleArray& dist;
|
||||
// BlobIDstruct last_id, new_index, new_id;
|
||||
// BlobIDList new_list;
|
||||
// MPI_Comm newcomm;
|
||||
//};
|
||||
//
|
||||
|
||||
|
||||
//
|
||||
//
|
||||
//class BlobIdentificationWorkItem2: public ThreadPool::WorkItem
|
||||
//{
|
||||
//public:
|
||||
// BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
|
||||
// std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
|
||||
// BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
|
||||
// timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
|
||||
// phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
|
||||
// {
|
||||
// MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
|
||||
// }
|
||||
// ~BlobIdentificationWorkItem2() { MPI_Comm_free(&newcomm); }
|
||||
// virtual void run() {
|
||||
// // Compute the global blob id and compare to the previous version
|
||||
// PROFILE_START("Identify blobs maps",1);
|
||||
// const IntArray& ids = new_index->second;
|
||||
// static int max_id = -1;
|
||||
// new_id->first = new_index->first;
|
||||
// new_id->second = new_index->second;
|
||||
// if ( last_id.get()!=NULL ) {
|
||||
// // Compute the timestep-timestep map
|
||||
// const IntArray& old_ids = last_id->second;
|
||||
// ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,newcomm);
|
||||
// // Renumber the current timestep's ids
|
||||
// getNewIDs(map,max_id,*new_list);
|
||||
// renumberIDs(*new_list,new_id->second);
|
||||
// writeIDMap(map,timestep,id_map_filename);
|
||||
// } else {
|
||||
// max_id = -1;
|
||||
// ID_map_struct map(new_id->first);
|
||||
// getNewIDs(map,max_id,*new_list);
|
||||
// writeIDMap(map,timestep,id_map_filename);
|
||||
// }
|
||||
// PROFILE_STOP("Identify blobs maps",1);
|
||||
// }
|
||||
// virtual bool has_result() const { return false; }
|
||||
//private:
|
||||
// BlobIdentificationWorkItem2();
|
||||
// int timestep;
|
||||
// int Nx, Ny, Nz;
|
||||
// const RankInfoStruct& rank_info;
|
||||
// std::shared_ptr<const DoubleArray> phase;
|
||||
// const DoubleArray& dist;
|
||||
// BlobIDstruct last_id, new_index, new_id;
|
||||
// BlobIDList new_list;
|
||||
// MPI_Comm newcomm;
|
||||
//};
|
||||
//
|
||||
|
||||
|
||||
// Helper class to write the vis file from a thread
|
||||
class WriteVisWorkItem: public ThreadPool::WorkItem
|
||||
{
|
||||
public:
|
||||
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
|
||||
TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
|
||||
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_)
|
||||
{
|
||||
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
|
||||
}
|
||||
~WriteVisWorkItem() { MPI_Comm_free(&newcomm); }
|
||||
virtual void run() {
|
||||
PROFILE_START("Save Vis",1);
|
||||
ASSERT(visData[0].vars[0]->name=="phase");
|
||||
ASSERT(visData[0].vars[1]->name=="Pressure");
|
||||
ASSERT(visData[0].vars[2]->name=="SignDist");
|
||||
ASSERT(visData[0].vars[3]->name=="BlobID");
|
||||
Array<double>& PhaseData = visData[0].vars[0]->data;
|
||||
Array<double>& PressData = visData[0].vars[1]->data;
|
||||
Array<double>& SignData = visData[0].vars[2]->data;
|
||||
Array<double>& BlobData = visData[0].vars[3]->data;
|
||||
fillData.copy(Averages.SDn,PhaseData);
|
||||
fillData.copy(Averages.Press,PressData);
|
||||
fillData.copy(Averages.SDs,SignData);
|
||||
fillData.copy(Averages.Label_NWP,BlobData);
|
||||
IO::writeData( timestep, visData, newcomm );
|
||||
PROFILE_STOP("Save Vis",1);
|
||||
};
|
||||
virtual bool has_result() const { return false; }
|
||||
private:
|
||||
WriteVisWorkItem();
|
||||
int timestep;
|
||||
std::vector<IO::MeshDataStruct>& visData;
|
||||
TwoPhase& Averages;
|
||||
fillHalo<double>& fillData;
|
||||
MPI_Comm newcomm;
|
||||
};
|
||||
|
||||
|
||||
|
||||
// Helper class to run the analysis from within a thread
|
||||
// Note: Averages will be modified after the constructor is called
|
||||
class AnalysisWorkItem: public ThreadPool::WorkItem
|
||||
{
|
||||
public:
|
||||
AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_,
|
||||
BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
|
||||
type(type_), timestep(timestep_), Averages(Averages_),
|
||||
blob_ids(ids), id_list(id_list_), beta(beta_) { }
|
||||
~AnalysisWorkItem() { }
|
||||
virtual void run() {
|
||||
Averages.NumberComponents_NWP = blob_ids->first;
|
||||
Averages.Label_NWP = blob_ids->second;
|
||||
Averages.Label_NWP_map = *id_list;
|
||||
Averages.NumberComponents_WP = 1;
|
||||
Averages.Label_WP.fill(0.0);
|
||||
if ( (type&CopyPhaseIndicator) != 0 ) {
|
||||
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
|
||||
}
|
||||
if ( (type&ComputeAverages) != 0 ) {
|
||||
PROFILE_START("Compute dist",1);
|
||||
Averages.Initialize();
|
||||
Averages.ComputeDelPhi();
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
|
||||
Averages.UpdateMeshValues();
|
||||
Averages.ComputeLocal();
|
||||
Averages.Reduce();
|
||||
Averages.PrintAll(timestep);
|
||||
Averages.Initialize();
|
||||
Averages.ComponentAverages();
|
||||
Averages.SortBlobs();
|
||||
Averages.PrintComponents(timestep);
|
||||
PROFILE_STOP("Compute dist",1);
|
||||
}
|
||||
}
|
||||
virtual bool has_result() const { return false; }
|
||||
private:
|
||||
AnalysisWorkItem();
|
||||
AnalysisType type;
|
||||
int timestep;
|
||||
TwoPhase& Averages;
|
||||
BlobIDstruct blob_ids;
|
||||
BlobIDList id_list;
|
||||
double beta;
|
||||
};
|
||||
|
||||
|
||||
// Function to start the analysis
|
||||
void run_analysis( int timestep, int restart_interval,
|
||||
const RankInfoStruct& rank_info, TwoPhase& Averages,
|
||||
BlobIDstruct& last_ids, BlobIDstruct& last_index, BlobIDList& last_id_map,
|
||||
int Nx, int Ny, int Nz, bool pBC, double err,
|
||||
double *Pressure, const double *Velocity,
|
||||
const char *ID, const double *f_even, const double *f_odd,
|
||||
const char *LocalRestartFile, std::vector<IO::MeshDataStruct>& visData, fillHalo<double>& fillData,
|
||||
ThreadPool& tpool, AnalysisWaitIdStruct& wait )
|
||||
{
|
||||
|
||||
int N = Nx*Ny*Nz;
|
||||
|
||||
/*
|
||||
* 2
|
||||
*/
|
||||
|
||||
// Determin the analysis we want to perform
|
||||
AnalysisType type = AnalyzeNone;
|
||||
// if ( timestep%ANALYSIS_INTERVAL + 5 == ANALYSIS_INTERVAL ) {
|
||||
// // Copy the phase indicator field for the earlier timestep
|
||||
// type = static_cast<AnalysisType>( type | CopyPhaseIndicator );
|
||||
// }
|
||||
// if ( timestep%BLOBID_INTERVAL == 0 ) {
|
||||
// // Identify blobs and update global ids in time
|
||||
// type = static_cast<AnalysisType>( type | IdentifyBlobs );
|
||||
// }
|
||||
//
|
||||
// if ( timestep%ANALYSIS_INTERVAL == 0 ) {
|
||||
// // Copy the averages to the CPU (and identify blobs)
|
||||
// type = static_cast<AnalysisType>( type | CopySimState );
|
||||
// type = static_cast<AnalysisType>( type | IdentifyBlobs );
|
||||
// }
|
||||
// if ( timestep%ANALYSIS_INTERVAL == 5 ) {
|
||||
// // Run the analysis
|
||||
// type = static_cast<AnalysisType>( type | ComputeAverages );
|
||||
// }
|
||||
// if (timestep%restart_interval == 0) {
|
||||
// // Write the restart file
|
||||
// type = static_cast<AnalysisType>( type | CreateRestart );
|
||||
// }
|
||||
if (timestep%restart_interval == 0) {
|
||||
// Write the visualization data
|
||||
type = static_cast<AnalysisType>( type | WriteVis );
|
||||
type = static_cast<AnalysisType>( type | CopySimState );
|
||||
type = static_cast<AnalysisType>( type | IdentifyBlobs );
|
||||
}
|
||||
|
||||
// Return if we are not doing anything
|
||||
if ( type == AnalyzeNone )
|
||||
return;
|
||||
|
||||
/*
|
||||
* 3
|
||||
*/
|
||||
|
||||
PROFILE_START("start_analysis");
|
||||
|
||||
// Copy the appropriate variables to the host (so we can spawn new threads)
|
||||
ScaLBL_DeviceBarrier();
|
||||
|
||||
/*
|
||||
* 4
|
||||
*/
|
||||
|
||||
PROFILE_START("Copy data to host",1);
|
||||
|
||||
std::shared_ptr<DoubleArray> phase;
|
||||
// if ( (type&CopyPhaseIndicator)!=0 || (type&ComputeAverages)!=0 || (type&CopySimState)!=0 || (type&IdentifyBlobs)!=0 )
|
||||
// {
|
||||
// phase = std::shared_ptr<DoubleArray>(new DoubleArray(Nx,Ny,Nz));
|
||||
// ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
|
||||
// }
|
||||
|
||||
|
||||
// if ( (type&CopyPhaseIndicator)!=0 ) {
|
||||
// memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
|
||||
// //Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
|
||||
// }
|
||||
|
||||
|
||||
// if ( (type&ComputeAverages)!=0 ) {
|
||||
// memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
|
||||
// //Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tminus);
|
||||
// }
|
||||
|
||||
|
||||
if ( (type&CopySimState) != 0 ) {
|
||||
// Copy the members of Averages to the cpu (phase was copied above)
|
||||
// Wait
|
||||
PROFILE_START("Copy-Pressure",1);
|
||||
ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
ScaLBL_DeviceBarrier();
|
||||
PROFILE_STOP("Copy-Pressure",1);
|
||||
PROFILE_START("Copy-Wait",1);
|
||||
tpool.wait(wait.analysis);
|
||||
tpool.wait(wait.vis); // Make sure we are done using analysis before modifying
|
||||
PROFILE_STOP("Copy-Wait",1);
|
||||
PROFILE_START("Copy-State",1);
|
||||
// memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
|
||||
PROFILE_STOP("Copy-State",1);
|
||||
}
|
||||
|
||||
std::shared_ptr<double> cDen, cDistEven, cDistOdd;
|
||||
if ( (type&CreateRestart) != 0 ) {
|
||||
// Copy restart data to the CPU
|
||||
cDen = std::shared_ptr<double>(new double[2*N],DeleteArray<double>);
|
||||
cDistEven = std::shared_ptr<double>(new double[10*N],DeleteArray<double>);
|
||||
cDistOdd = std::shared_ptr<double>(new double[9*N],DeleteArray<double>);
|
||||
ScaLBL_CopyToHost(cDistEven.get(),f_even,10*N*sizeof(double));
|
||||
ScaLBL_CopyToHost(cDistOdd.get(),f_odd,9*N*sizeof(double));
|
||||
// ScaLBL_CopyToHost(cDen.get(),Den,2*N*sizeof(double));
|
||||
}
|
||||
|
||||
|
||||
PROFILE_STOP("Copy data to host",1);
|
||||
|
||||
|
||||
/*
|
||||
* 5
|
||||
*/
|
||||
|
||||
//
|
||||
// // Spawn threads to do blob identification work
|
||||
// if ( (type&IdentifyBlobs)!=0 ) {
|
||||
// BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
|
||||
// BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
|
||||
// BlobIDList new_list(new std::vector<BlobIDType>());
|
||||
// ThreadPool::WorkItem *work1 = new BlobIdentificationWorkItem1(timestep,
|
||||
// Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list);
|
||||
// ThreadPool::WorkItem *work2 = new BlobIdentificationWorkItem2(timestep,
|
||||
// Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list);
|
||||
// work1->add_dependency(wait.blobID);
|
||||
// work2->add_dependency(tpool.add_work(work1));
|
||||
// wait.blobID = tpool.add_work(work2);
|
||||
// last_index = new_index;
|
||||
// last_ids = new_ids;
|
||||
// last_id_map = new_list;
|
||||
// }
|
||||
//
|
||||
|
||||
/*
|
||||
* 6
|
||||
*/
|
||||
|
||||
// // Spawn threads to do the analysis work
|
||||
// if ( (type&ComputeAverages) != 0 ) {
|
||||
// double beta = 0.0;
|
||||
// ThreadPool::WorkItem *work = new AnalysisWorkItem(type,timestep,Averages,last_index,last_id_map,beta);
|
||||
// work->add_dependency(wait.blobID);
|
||||
// work->add_dependency(wait.analysis);
|
||||
// work->add_dependency(wait.vis); // Make sure we are done using analysis before modifying
|
||||
// wait.analysis = tpool.add_work(work);
|
||||
// }
|
||||
//
|
||||
|
||||
/*
|
||||
* 7
|
||||
*/
|
||||
|
||||
// Spawn a thread to write the restart file
|
||||
if ( (type&CreateRestart) != 0 ) {
|
||||
int rank = MPI_WORLD_RANK();
|
||||
|
||||
// Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
|
||||
tpool.wait(wait.restart);
|
||||
// Retain the timestep associated with the restart files
|
||||
if (rank==0){
|
||||
FILE *Rst = fopen("Restart.txt","w");
|
||||
fprintf(Rst,"%i\n",timestep+5);
|
||||
fclose(Rst);
|
||||
}
|
||||
// Write the restart file (using a seperate thread)
|
||||
WriteRestartWorkItem *work = new WriteRestartWorkItem(LocalRestartFile,cDen,cDistEven,cDistOdd,N);
|
||||
work->add_dependency(wait.restart);
|
||||
wait.restart = tpool.add_work(work);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* 8
|
||||
*/
|
||||
|
||||
// Save the results for visualization
|
||||
if ( (type&CreateRestart) != 0 ) {
|
||||
// Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
|
||||
tpool.wait(wait.vis);
|
||||
// Write the vis files
|
||||
// ThreadPool::WorkItem *work = new WriteVisWorkItem( timestep, visData, Averages, fillData );
|
||||
WriteVisWorkItem *work = new WriteVisWorkItem( timestep, visData, Averages, fillData );
|
||||
work->add_dependency(wait.blobID);
|
||||
work->add_dependency(wait.analysis);
|
||||
work->add_dependency(wait.vis);
|
||||
wait.vis = tpool.add_work(work);
|
||||
}
|
||||
|
||||
|
||||
PROFILE_STOP("start_analysis");
|
||||
printf("Returning from the run_analysis function...\n");
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
//if (pBC) {
|
||||
//err = fabs(sat_w - sat_w_previous);
|
||||
//sat_w_previous = sat_w;
|
||||
//if (rank==0){
|
||||
// printf("Timestep %i: change in saturation since last checkpoint is %f \n",timestep,err);
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/* #ifdef USE_CUDA
|
||||
if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
|
||||
// Keep a few blob identifications queued up to keep the processors busy,
|
||||
// allowing us to track the blobs as fast as possible
|
||||
// Add more detailed estimates of the update frequency required to track blobs
|
||||
type = static_cast<AnalysisType>( type | IdentifyBlobs );
|
||||
}
|
||||
#endif
|
||||
*/
|
||||
Reference in New Issue
Block a user