diff --git a/analysis/FreeEnergy.cpp b/analysis/FreeEnergy.cpp new file mode 100644 index 00000000..6a641a95 --- /dev/null +++ b/analysis/FreeEnergy.cpp @@ -0,0 +1,181 @@ +#include "analysis/FreeEnergy.h" + +FreeEnergyAnalyzer::FreeEnergyAnalyzer(std::shared_ptr dm): + Dm(dm) +{ + + Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz; + Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0; + + ChemicalPotential.resize(Nx,Ny,Nz); ChemicalPotential.fill(0); + Phi.resize(Nx,Ny,Nz); Phi.fill(0); + Pressure.resize(Nx,Ny,Nz); Pressure.fill(0); + Rho.resize(Nx,Ny,Nz); Rho.fill(0); + Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field + Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0); + Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0); + SDs.resize(Nx,Ny,Nz); SDs.fill(0); + + if (Dm->rank()==0){ + bool WriteHeader=false; + TIMELOG = fopen("free.csv","r"); + if (TIMELOG != NULL) + fclose(TIMELOG); + else + WriteHeader=true; + + TIMELOG = fopen("free.csv","a+"); + if (WriteHeader) + { + // If timelog is empty, write a short header to list the averages + //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); + fprintf(TIMELOG,"timestep\n"); + } + } + +} + +FreeEnergyAnalyzer::~FreeEnergyAnalyzer(){ + if (Dm->rank()==0){ + fclose(TIMELOG); + } +} + +void FreeEnergyAnalyzer::SetParams(){ + +} + +void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){ + + int i,j,k; + + if (Dm->rank()==0){ + fprintf(TIMELOG,"%i ",timestep); + /*for (int ion=0; ion input_db, int timestep){ + + auto vis_db = input_db->getDatabase( "Visualization" ); + char VisName[40]; + + std::vector visData; + fillHalo fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1); + + IO::initialize("","silo","false"); + // Create the MeshDataStruct + visData.resize(1); + + visData[0].meshName = "domain"; + visData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); + auto VisPhase = std::make_shared(); + auto VisPressure = std::make_shared(); + auto VisChemicalPotential = std::make_shared(); + auto VxVar = std::make_shared(); + auto VyVar = std::make_shared(); + auto VzVar = std::make_shared(); + + + if (vis_db->getWithDefault( "save_phase_field", true )){ + VisPhase->name = "Phase"; + VisPhase->type = IO::VariableType::VolumeVariable; + VisPhase->dim = 1; + VisPhase->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VisPhase); + } + + if (vis_db->getWithDefault( "save_potential", true )){ + + VisPressure->name = "Pressure"; + VisPressure->type = IO::VariableType::VolumeVariable; + VisPressure->dim = 1; + VisPressure->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VisPressure); + + VisChemicalPotential->name = "ChemicalPotential"; + VisChemicalPotential->type = IO::VariableType::VolumeVariable; + VisChemicalPotential->dim = 1; + VisChemicalPotential->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VisChemicalPotential); + } + + if (vis_db->getWithDefault( "save_velocity", false )){ + VxVar->name = "Velocity_x"; + VxVar->type = IO::VariableType::VolumeVariable; + VxVar->dim = 1; + VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VxVar); + VyVar->name = "Velocity_y"; + VyVar->type = IO::VariableType::VolumeVariable; + VyVar->dim = 1; + VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VyVar); + VzVar->name = "Velocity_z"; + VzVar->type = IO::VariableType::VolumeVariable; + VzVar->dim = 1; + VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VzVar); + } + + if (vis_db->getWithDefault( "save_phase", true )){ + ASSERT(visData[0].vars[0]->name=="Phase"); + LeeModel.getPhase(Phi); + Array& PhaseData = visData[0].vars[0]->data; + fillData.copy(Phi,PhaseData); + } + + if (vis_db->getWithDefault( "save_potential", true )){ + ASSERT(visData[0].vars[1]->name=="Pressure"); + LeeModel.getPotential(Pressure, ChemicalPotential); + Array& PressureData = visData[0].vars[1]->data; + fillData.copy(Pressure,PressureData); + + ASSERT(visData[0].vars[2]->name=="ChemicalPotential"); + Array& ChemicalPotentialData = visData[0].vars[2]->data; + fillData.copy(ChemicalPotential,ChemicalPotentialData); + } + + if (vis_db->getWithDefault( "save_velocity", false )){ + ASSERT(visData[0].vars[3]->name=="Velocity_x"); + ASSERT(visData[0].vars[4]->name=="Velocity_y"); + ASSERT(visData[0].vars[5]->name=="Velocity_z"); + LeeModel.getVelocity(Vel_x,Vel_y,Vel_z); + Array& VelxData = visData[0].vars[3]->data; + Array& VelyData = visData[0].vars[4]->data; + Array& VelzData = visData[0].vars[5]->data; + fillData.copy(Vel_x,VelxData); + fillData.copy(Vel_y,VelyData); + fillData.copy(Vel_z,VelzData); + } + + if (vis_db->getWithDefault( "write_silo", true )) + IO::writeData( timestep, visData, Dm->Comm ); + +/* if (vis_db->getWithDefault( "save_8bit_raw", true )){ + char CurrentIDFilename[40]; + sprintf(CurrentIDFilename,"id_t%d.raw",timestep); + Averages.AggregateLabels(CurrentIDFilename); + } +*/ +} diff --git a/analysis/FreeEnergy.h b/analysis/FreeEnergy.h new file mode 100644 index 00000000..fbb1ba31 --- /dev/null +++ b/analysis/FreeEnergy.h @@ -0,0 +1,54 @@ +/* + * averaging tools for electrochemistry + */ + +#ifndef FreeEnergyAnalyzer_INC +#define FreeEnergyAnalyzer_INC + +#include +#include "common/Domain.h" +#include "common/Utilities.h" +#include "common/MPI.h" +#include "common/Communication.h" +#include "analysis/analysis.h" +#include "analysis/distance.h" +#include "analysis/Minkowski.h" +#include "analysis/SubPhase.h" +#include "IO/MeshDatabase.h" +#include "IO/Reader.h" +#include "IO/Writer.h" +#include "models/FreeLeeModel.h" + +class FreeEnergyAnalyzer{ +public: + std::shared_ptr Dm; + double Volume; + // input variables + double rho_n, rho_w; + double nu_n, nu_w; + double gamma_wn, beta; + double Fx, Fy, Fz; + + //........................................................................... + int Nx,Ny,Nz; + DoubleArray Rho; + DoubleArray Phi; + DoubleArray ChemicalPotential; + DoubleArray Pressure; + DoubleArray Vel_x; + DoubleArray Vel_y; + DoubleArray Vel_z; + DoubleArray SDs; + + FreeEnergyAnalyzer(std::shared_ptr Dm); + ~FreeEnergyAnalyzer(); + + void SetParams(); + void Basic( ScaLBL_FreeLeeModel &LeeModel, int timestep); + void WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr input_db, int timestep); + +private: + FILE *TIMELOG; +}; +#endif + diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 37f58d0c..f21767dd 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -702,12 +702,14 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array &id, if (rank == 0) printf(" delta=%f, growth=%f, max. displacement = %f \n",morph_delta, GrowthEstimate, MAX_DISPLACEMENT); // Now adjust morph_delta - double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious); - GrowthPrevious = GrowthEstimate; - morph_delta_previous = morph_delta; - morph_delta += step_size; + if (fabs(GrowthEstimate - GrowthPrevious) > 0.0) { + double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious); + GrowthPrevious = GrowthEstimate; + morph_delta_previous = morph_delta; + morph_delta += step_size; + } if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0; - + //MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25); if (morph_delta > 0.0 ){ diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index f43a26ff..ab40ae4c 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -706,6 +706,139 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStru } + // Initialize the comms + for ( int i = 0; i < 1024; i++ ) + d_comm_used[i] = false; + // Initialize the threads + int N_threads = db->getWithDefault( "N_threads", 4 ); + auto method = db->getWithDefault( "load_balance", "default" ); + createThreads( method, N_threads ); +} + +runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel) +/* std::shared_ptr input_db, const RankInfoStruct &rank_info, + std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, int Np, + bool Regular, IntArray Map ) + : d_Np( Np ), + d_regular( Regular ), + d_rank_info( rank_info ), + d_Map( Map ), + d_comm( Dm->Comm.dup() ), + d_ScaLBL_Comm( ScaLBL_Comm )*/ +{ + + d_comm = ColorModel.Dm->Comm.dup(); + d_Np = ColorModel.Np; + bool Regular = false; + + auto input_db = ColorModel.db; + auto db = input_db->getDatabase( "Analysis" ); + auto vis_db = input_db->getDatabase( "Visualization" ); + + // Ids of work items to use for dependencies + ThreadPool::thread_id_t d_wait_blobID; + ThreadPool::thread_id_t d_wait_analysis; + ThreadPool::thread_id_t d_wait_vis; + ThreadPool::thread_id_t d_wait_restart; + ThreadPool::thread_id_t d_wait_subphase; + + char rankString[20]; + sprintf( rankString, "%05d", ColorModel.Dm->rank() ); + d_n[0] = ColorModel.Dm->Nx - 2; + d_n[1] = ColorModel.Dm->Ny - 2; + d_n[2] = ColorModel.Dm->Nz - 2; + d_N[0] = ColorModel.Dm->Nx; + d_N[1] = ColorModel.Dm->Ny; + d_N[2] = ColorModel.Dm->Nz; + + d_restart_interval = db->getScalar( "restart_interval" ); + d_analysis_interval = db->getScalar( "analysis_interval" ); + d_subphase_analysis_interval = INT_MAX; + d_visualization_interval = INT_MAX; + d_blobid_interval = INT_MAX; + if ( db->keyExists( "blobid_interval" ) ) { + d_blobid_interval = db->getScalar( "blobid_interval" ); + } + if ( db->keyExists( "visualization_interval" ) ) { + d_visualization_interval = db->getScalar( "visualization_interval" ); + } + if ( db->keyExists( "subphase_analysis_interval" ) ) { + d_subphase_analysis_interval = db->getScalar( "subphase_analysis_interval" ); + } + + auto restart_file = db->getScalar( "restart_file" ); + d_restartFile = restart_file + "." + rankString; + + + d_rank = d_comm.getRank(); + writeIDMap( ID_map_struct(), 0, id_map_filename ); + // Initialize IO for silo + IO::initialize( "", "silo", "false" ); + // Create the MeshDataStruct + d_meshData.resize( 1 ); + + d_meshData[0].meshName = "domain"; + d_meshData[0].mesh = std::make_shared( + d_rank_info, d_n[0], d_n[1], d_n[2], ColorModel.Dm->Lx, ColorModel.Dm->Ly, ColorModel.Dm->Lz ); + auto PhaseVar = std::make_shared(); + auto PressVar = std::make_shared(); + auto VxVar = std::make_shared(); + auto VyVar = std::make_shared(); + auto VzVar = std::make_shared(); + auto SignDistVar = std::make_shared(); + auto BlobIDVar = std::make_shared(); + + if ( vis_db->getWithDefault( "save_phase_field", true ) ) { + PhaseVar->name = "phase"; + PhaseVar->type = IO::VariableType::VolumeVariable; + PhaseVar->dim = 1; + PhaseVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( PhaseVar ); + } + + if ( vis_db->getWithDefault( "save_pressure", false ) ) { + PressVar->name = "Pressure"; + PressVar->type = IO::VariableType::VolumeVariable; + PressVar->dim = 1; + PressVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( PressVar ); + } + + if ( vis_db->getWithDefault( "save_velocity", false ) ) { + VxVar->name = "Velocity_x"; + VxVar->type = IO::VariableType::VolumeVariable; + VxVar->dim = 1; + VxVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VxVar ); + VyVar->name = "Velocity_y"; + VyVar->type = IO::VariableType::VolumeVariable; + VyVar->dim = 1; + VyVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VyVar ); + VzVar->name = "Velocity_z"; + VzVar->type = IO::VariableType::VolumeVariable; + VzVar->dim = 1; + VzVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VzVar ); + } + + if ( vis_db->getWithDefault( "save_distance", false ) ) { + SignDistVar->name = "SignDist"; + SignDistVar->type = IO::VariableType::VolumeVariable; + SignDistVar->dim = 1; + SignDistVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( SignDistVar ); + } + + if ( vis_db->getWithDefault( "save_connected_components", false ) ) { + BlobIDVar->name = "BlobID"; + BlobIDVar->type = IO::VariableType::VolumeVariable; + BlobIDVar->dim = 1; + BlobIDVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( BlobIDVar ); + } + + // Initialize the comms for ( int i = 0; i < 1024; i++ ) d_comm_used[i] = false; diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index a82c4ba0..c7c4ce71 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -7,6 +7,7 @@ #include "common/Communication.h" #include "common/ScaLBL.h" #include "threadpool/thread_pool.h" +#include "models/ColorModel.h" #include @@ -31,6 +32,8 @@ public: runAnalysis( std::shared_ptr db, const RankInfoStruct &rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr dm, int Np, bool Regular, IntArray Map ); + + runAnalysis( ScaLBL_ColorModel &ColorModel); //! Destructor ~runAnalysis(); diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index dcadb08e..4726bae6 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -516,9 +516,9 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis n = k*Nx*Ny+j*Nx+i; if (id[n] > 0){ // Counts for the six faces - if (i>0 && i<=width) Map(n)=idx++; - else if (j>0 && j<=width) Map(n)=idx++; - else if (k>0 && k<=width) Map(n)=idx++; + if (i>0 && i<=width) Map(n)=idx++; + else if (j>0 && j<=width) Map(n)=idx++; + else if (k>0 && k<=width) Map(n)=idx++; else if (i>Nx-width-2 && iNy-width-2 && jNz-width-2 && k0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + // q = 1 + //nread = neighborList[n+Np]; + Aq[nr2] = a1; + Bq[nr2] = b1; + // q=2 + //nread = neighborList[n]; + Aq[nr1] = a2; + Bq[nr1] = b2; + + //............................................... + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + // q = 3 + //nread = neighborList[n+3*Np]; + Aq[nr4] = a1; + Bq[nr4] = b1; + // q = 4 + //nread = neighborList[n+2*Np]; + Aq[nr3] = a2; + Bq[nr3] = b2; + + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + // q = 5 + //nread = neighborList[n+5*Np]; + Aq[nr6] = a1; + Bq[nr6] = b1; + // q = 6 + //nread = neighborList[n+4*Np]; + Aq[nr5] = a2; + Bq[nr5] = b2; + //............................................... + } +} + +extern "C" void ScaLBL_D3Q7_AAeven_Color(int *Map, double *Aq, double *Bq, double *Den, + double *Phi, double *ColorGrad, double *Vel, double rhoA, double rhoB, double beta, int start, int finish, int Np){ + + double nA,nB; // number density + double a1,b1,a2,b2,nAB,delta; + double C,nx,ny,nz; //color gradient magnitude and direction + double ux,uy,uz; + double phi; + // Instantiate mass transport distributions + // Stationary value - distribution 0 + for (int n=start; n0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta; + + Aq[1*Np+n] = a1; + Bq[1*Np+n] = b1; + Aq[2*Np+n] = a2; + Bq[2*Np+n] = b2; + + //............................................... + // q = 2 + // Cq = {0,1,0} + delta = beta*nA*nB*nAB*0.1111111111111111*ny; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta; + + Aq[3*Np+n] = a1; + Bq[3*Np+n] = b1; + Aq[4*Np+n] = a2; + Bq[4*Np+n] = b2; + //............................................... + // q = 4 + // Cq = {0,0,1} + delta = beta*nA*nB*nAB*0.1111111111111111*nz; + if (!(nA*nB*nAB>0)) delta=0; + a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; + b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; + a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta; + b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta; + + Aq[5*Np+n] = a1; + Bq[5*Np+n] = b1; + Aq[6*Np+n] = a2; + Bq[6*Np+n] = b2; + //............................................... + + } +} + + extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np){ - int idx,n,nread; + int idx,nread; double fq,nA,nB; for (int n=start; n +#include #define STOKES @@ -70,6 +71,8 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double gqbar[17*Np+n] = 0.0277777777777778*(p-0.5*(Fy-Fz)); ; //double(100*n)+17.f; gqbar[18*Np+n] = 0.0277777777777778*(p-0.5*(-Fy+Fz));; //double(100*n)+18.f; } + + } extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, @@ -101,7 +104,8 @@ extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, doubl nz = nz/ColorMag_temp; theta = M*cs2_inv*(1-4.0*phi*phi)/W; - + theta = 0; // try more diffusive initial condition + hq[0*Np+idx]=0.3333333333333333*(phi); hq[1*Np+idx]=0.1111111111111111*(phi+theta*nx); hq[2*Np+idx]=0.1111111111111111*(phi-theta*nx); @@ -116,7 +120,7 @@ extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, doubl extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){ - int idx,n,nread; + int idx,nread; double fq,phi; for (int n=start; n +#include +#include #define STOKES @@ -185,10 +187,16 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, } } -__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ - int idx,n; - double fq,phi; +__global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ + + int n,idx,nr1,nr2,nr3,nr4,nr5,nr6; + double h0,h1,h2,h3,h4,h5,h6; + double nx,ny,nz,C; + double ux,uy,uz; + double phi; + double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 + double factor = 1.0; // for (int n=start; n>>( gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init: %s \n",cudaGetErrorString(err)); + } } extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){ - + + dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<>>( gqbar, Fx, Fy, Fz, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init: %s \n",cudaGetErrorString(err)); + } } extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ -} -extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ + dvc_ScaLBL_FreeLeeModel_PhaseField_Init<<>>(Map, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_FreeLeeModel_PhaseField_Init: %s \n",cudaGetErrorString(err)); + } + } - -extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ - +extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np) +{ + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField<<>>(neighborList, Map, hq, Den, Phi, ColorGrad, Vel, + rhoA, rhoB, tauM, W, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAodd_FreeLee_PhaseField: %s \n",cudaGetErrorString(err)); + } } -extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, - double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, +extern "C" void ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ + + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField<<>>( Map, hq, Den, Phi, ColorGrad, Vel, rhoA, rhoB, tauM, W, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAeven_FreeLee_PhaseField: %s \n",cudaGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){ + + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q7_ComputePhaseField, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q7_ComputePhaseField<<>>( Map, hq, Den, Phi, rhoA, rhoB, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_ComputePhaseField: %s \n",cudaGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, + double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -} - + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel<<>>(neighborList, Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad, + rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel: %s \n",cudaGetErrorString(err)); + } +} -extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, - double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, +extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, + double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel<<>>(Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad, + rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel: %s \n",cudaGetErrorString(err)); + } + } extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ - + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK<<>>(neighborList, dist, Vel, Pressure, + tau, rho0, Fx, Fy, Fz, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK: %s \n",cudaGetErrorString(err)); + } } extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ - + + cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK, cudaFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK<<>>(dist, Vel, Pressure, + tau, rho0, Fx, Fy, Fz, start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK: %s \n",cudaGetErrorString(err)); + } } + +extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np){ +} \ No newline at end of file diff --git a/hip/FreeLee.cu b/hip/FreeLee.cu index 558bd2f1..09bc8689 100644 --- a/hip/FreeLee.cu +++ b/hip/FreeLee.cu @@ -1,11 +1,12 @@ #include +#include #include "hip/hip_runtime.h" +#define STOKES + #define NBLOCKS 1024 #define NTHREADS 256 -#define STOKES - __global__ void dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np) { int n; @@ -186,10 +187,16 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, } } -__global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ - int idx,n; - double fq,phi; +__global__ void dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ + + int n,idx,nr1,nr2,nr3,nr4,nr5,nr6; + double h0,h1,h2,h3,h4,h5,h6; + double nx,ny,nz,C; + double ux,uy,uz; + double phi; + double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 + double factor = 1.0; // for (int n=start; n>>( gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init: %s \n",hipGetErrorString(err)); + } } extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){ - + + dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<>>( gqbar, Fx, Fy, Fz, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init: %s \n",hipGetErrorString(err)); + } } extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ -} -extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ + dvc_ScaLBL_FreeLeeModel_PhaseField_Init<<>>(Map, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_FreeLeeModel_PhaseField_Init: %s \n",hipGetErrorString(err)); + } + } - -extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, - double rhoA, double rhoB, int start, int finish, int Np){ - +extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np) +{ + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField<<>>(neighborList, Map, hq, Den, Phi, ColorGrad, Vel, + rhoA, rhoB, tauM, W, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAodd_FreeLee_PhaseField: %s \n",hipGetErrorString(err)); + } } -extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, - double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, +extern "C" void ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel, + double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField<<>>( Map, hq, Den, Phi, ColorGrad, Vel, rhoA, rhoB, tauM, W, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_AAeven_FreeLee_PhaseField: %s \n",hipGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){ + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_ComputePhaseField, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q7_ComputePhaseField<<>>( Map, hq, Den, Phi, rhoA, rhoB, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q7_ComputePhaseField: %s \n",hipGetErrorString(err)); + } +} + + +extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, + double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -} - + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel<<>>(neighborList, Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad, + rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel: %s \n",hipGetErrorString(err)); + } +} -extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, - double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, +extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, + double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel<<>>(Map, dist, Den, Phi, mu_phi, Vel, Pressure, ColorGrad, + rhoA, rhoB, tauA, tauB, kappa, beta, W, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel: %s \n",hipGetErrorString(err)); + } + } extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ - + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK<<>>(neighborList, dist, Vel, Pressure, + tau, rho0, Fx, Fy, Fz, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK: %s \n",hipGetErrorString(err)); + } } extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ - + + hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK, hipFuncCachePreferL1); + dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK<<>>(dist, Vel, Pressure, + tau, rho0, Fx, Fy, Fz, start, finish, Np); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK: %s \n",hipGetErrorString(err)); + } } + +extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np){ +} \ No newline at end of file diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 1fe00824..204fd1d6 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -531,6 +531,121 @@ void ScaLBL_ColorModel::Initialize(){ ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); } +double ScaLBL_ColorModel::Run(int returntime){ + int nprocs=nprocx*nprocy*nprocz; + + //************ MAIN ITERATION LOOP ***************************************/ + comm.barrier(); + PROFILE_START("Loop"); + //std::shared_ptr analysis_db; + bool Regular = false; + auto current_db = db->cloneDatabase(); + auto t1 = std::chrono::system_clock::now(); + int START_TIMESTEP = timestep; + int EXIT_TIMESTEP = min(timestepMax,returntime); + while (timestep < EXIT_TIMESTEP ) { + //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } + PROFILE_START("Update"); + // *************ODD TIMESTEP************* + timestep++; + // Compute the Phase indicator field + // Read for Aq, Bq happens in this routine (requires communication) + ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + if (BoundaryCondition > 0 && BoundaryCondition < 5){ + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + // Halo exchange for phase field + ScaLBL_Comm_Regular->SendHalo(Phi); + + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set BCs + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + + // *************EVEN TIMESTEP************* + timestep++; + // Compute the Phase indicator field + ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + // Halo exchange for phase field + if (BoundaryCondition > 0 && BoundaryCondition < 5){ + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + ScaLBL_Comm_Regular->SendHalo(Phi); + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, + alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + //************************************************************************ + } + PROFILE_STOP("Update"); + + PROFILE_STOP("Loop"); + PROFILE_SAVE("lbpm_color_simulator",1); + //************************************************************************ + // Compute the walltime per timestep + auto t2 = std::chrono::system_clock::now(); + double cputime = std::chrono::duration( t2 - t1 ).count() / (timestep - START_TIMESTEP); + // Performance obtained from each node + double MLUPS = double(Np)/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); + return(MLUPS); + MLUPS *= nprocs; + +} + void ScaLBL_ColorModel::Run(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); @@ -580,7 +695,6 @@ void ScaLBL_ColorModel::Run(){ if (color_db->keyExists( "krA_morph_factor" )){ KRA_MORPH_FACTOR = color_db->getScalar( "krA_morph_factor" ); } - /* defaults for simulation protocols */ auto protocol = color_db->getWithDefault( "protocol", "none" ); if (protocol == "image sequence"){ @@ -625,7 +739,7 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "seed_water" )){ seed_water = analysis_db->getScalar( "seed_water" ); if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water); - USE_SEED = true; + ASSERT(protocol == "seed water"); } if (analysis_db->keyExists( "morph_delta" )){ morph_delta = analysis_db->getScalar( "morph_delta" ); @@ -656,7 +770,6 @@ void ScaLBL_ColorModel::Run(){ MAX_MORPH_TIMESTEPS = analysis_db->getScalar( "max_morph_timesteps" ); } - if (rank==0){ printf("********************************************************\n"); if (protocol == "image sequence"){ @@ -1320,7 +1433,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta double vF = 0.f; double vS = 0.f; double delta_volume; - double WallFactor = 0.0; + double WallFactor = 1.0; bool USE_CONNECTED_NWP = false; DoubleArray phase(Nx,Ny,Nz); @@ -1343,6 +1456,11 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta } } double volume_initial = Dm->Comm.sumReduce( count); + double PoreVolume = Dm->Volume*Dm->Porosity(); + /*ensure target isn't an absurdly small fraction of pore volume */ + if (volume_initial < target_delta_volume*PoreVolume){ + volume_initial = target_delta_volume*PoreVolume; + } /* sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank); FILE *INPUT = fopen(LocalRankFilename,"wb"); diff --git a/models/ColorModel.h b/models/ColorModel.h index b2a9c1d1..7d3c858a 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -16,6 +16,10 @@ Implementation of color lattice boltzmann model #include "ProfilerApp.h" #include "threadpool/thread_pool.h" + +#ifndef ScaLBL_ColorModel_INC +#define ScaLBL_ColorModel_INC + class ScaLBL_ColorModel{ public: ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM); @@ -29,6 +33,7 @@ public: void Create(); void Initialize(); void Run(); + double Run(int returntime); void WriteDebug(); void getPhaseField(DoubleArray &f); @@ -99,4 +104,5 @@ private: int timestep; int timestep_previous; }; +#endif diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index 5a6805b6..428db40f 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -10,9 +10,9 @@ color lattice boltzmann model #include ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM): -rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauM(0),rhoA(0),rhoB(0),W(0),gamma(0),kappa(0),beta(0), +rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(2),tauA(1.0),tauB(1.0),tauM(1.0),rhoA(1.0),rhoB(1.0),W(5.0),gamma(0.001),kappa(0.0075),beta(0.0024), Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0), -tau(0),rho0(0), +tau(1.0),rho0(1.0), Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) { @@ -20,6 +20,45 @@ Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0), ScaLBL_FreeLeeModel::~ScaLBL_FreeLeeModel(){ } + + +void ScaLBL_FreeLeeModel::getPhase(DoubleArray &PhaseValues){ + + DoubleArray PhaseWideHalo(Nxh,Nyh,Nzh); + ScaLBL_CopyToHost(PhaseWideHalo.data(), Phi, sizeof(double)*Nh); + + // use halo width = 1 for analysis data + for (int k=1; kRegularLayout(Map,Pressure,PressureValues); + ScaLBL_Comm->Barrier(); comm.barrier(); + + ScaLBL_Comm->RegularLayout(Map,mu_phi,MuValues); + ScaLBL_Comm->Barrier(); comm.barrier(); + +} + +void ScaLBL_FreeLeeModel::getVelocity(DoubleArray &Vel_x, DoubleArray &Vel_y, DoubleArray &Vel_z){ + + ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Vel_x); + ScaLBL_Comm->Barrier(); comm.barrier(); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Vel_y); + ScaLBL_Comm->Barrier(); comm.barrier(); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Vel_z); + ScaLBL_Comm->Barrier(); comm.barrier(); +} + void ScaLBL_FreeLeeModel::ReadParams(string filename){ // read the input database db = std::make_shared( filename ); @@ -529,6 +568,34 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad() //fwrite(phase,8,Nh,OUTFILE); //fclose(OUTFILE); + DoubleArray PhaseField(Nx,Ny,Nz); + FILE *OUTFILE; + ScaLBL_Comm->RegularLayout(Map,mu_phi_host,PhaseField); + sprintf(LocalRankFilename,"Chem_Init.%05i.raw",rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,OUTFILE); + fclose(OUTFILE); + + ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[0],PhaseField); + FILE *CGX_FILE; + sprintf(LocalRankFilename,"Gradient_X_Init.%05i.raw",rank); + CGX_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGX_FILE); + fclose(CGX_FILE); + + ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[Np],PhaseField); + FILE *CGY_FILE; + sprintf(LocalRankFilename,"Gradient_Y_Init.%05i.raw",rank); + CGY_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGY_FILE); + fclose(CGY_FILE); + + ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[2*Np],PhaseField); + FILE *CGZ_FILE; + sprintf(LocalRankFilename,"Gradient_Z_Init.%05i.raw",rank); + CGZ_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,CGZ_FILE); + fclose(CGZ_FILE); delete [] phase; delete [] ColorGrad_host; @@ -709,21 +776,17 @@ void ScaLBL_FreeLeeModel::Initialize_SingleFluid(){ } } -void ScaLBL_FreeLeeModel::Run_TwoFluid(){ +double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){ int nprocs=nprocx*nprocy*nprocz; - const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - if (rank==0){ - printf("********************************************************\n"); - printf("No. of timesteps: %i \n", timestepMax); - fflush(stdout); - } - + int START_TIME = timestep; + int EXIT_TIME = min(returntime, timestepMax); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); auto t1 = std::chrono::system_clock::now(); PROFILE_START("Loop"); - while (timestep < timestepMax ) { + + while (timestep < EXIT_TIME ) { //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } PROFILE_START("Update"); // *************ODD TIMESTEP************* @@ -732,24 +795,27 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){ // Compute the Phase indicator field // Read for hq happens in this routine (requires communication) ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np); // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL + // Halo exchange for phase field + ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_WideHalo->Send(Phi); + ScaLBL_Comm_WideHalo->Recv(Phi); if (BoundaryCondition > 0 && BoundaryCondition < 5){ //TODO to be revised + // Need to add BC for hq!!! ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } - // Halo exchange for phase field - ScaLBL_Comm_WideHalo->Send(Phi); + + ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL + ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, + kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM, - kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_WideHalo->Recv(Phi); ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); // Set BCs @@ -765,30 +831,34 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){ ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar); ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar); } - ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM, - kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); + ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, + kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + + // *************EVEN TIMESTEP************* timestep++; // Compute the Phase indicator field - ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL - ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMA + ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np); // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL // Halo exchange for phase field + ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_WideHalo->Send(Phi); + ScaLBL_Comm_WideHalo->Recv(Phi); if (BoundaryCondition > 0 && BoundaryCondition < 5){ ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } - ScaLBL_Comm_WideHalo->Send(Phi); - ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM, - kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_WideHalo->Recv(Phi); + ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL + + ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, + kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE ScaLBL_Comm->Barrier(); // Set boundary conditions @@ -804,8 +874,8 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){ ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar); ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar); } - ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM, - kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, + kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //************************************************************************ PROFILE_STOP("Update"); @@ -816,18 +886,11 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){ if (rank==0) printf("-------------------------------------------------------------------\n"); // Compute the walltime per timestep auto t2 = std::chrono::system_clock::now(); - double cputime = std::chrono::duration( t2 - t1 ).count() / timestep; + double cputime = std::chrono::duration( t2 - t1 ).count() / (EXIT_TIME-START_TIME); // Performance obtained from each node double MLUPS = double(Np)/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"); - - // ************************************************************************ + return MLUPS; } void ScaLBL_FreeLeeModel::Run_SingleFluid(){ @@ -878,6 +941,7 @@ void ScaLBL_FreeLeeModel::Run_SingleFluid(){ 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); + // *************EVEN TIMESTEP************* timestep++; //------------------------------------------------------------------------------------------------------------------- @@ -932,6 +996,32 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){ DoubleArray PhaseData(Nxh,Nyh,Nzh); //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); ScaLBL_CopyToHost(PhaseData.data(), Phi, sizeof(double)*Nh); + /* + IntArray MapData(Np); + ScaLBL_CopyToHost(MapData.data(), dvcMap, sizeof(int)*Np); + FILE *MAP; + sprintf(LocalRankFilename,"Map.%05i.raw",rank); + MAP = fopen(LocalRankFilename,"wb"); + fwrite(MapData.data(),4,Np,MAP); + fclose(MAP); + + FILE *NB; + //IntArray Neighbors(18,Np); + //ScaLBL_CopyToHost(Neighbors.data(), NeighborList, sizeof(int)*Np*18); + sprintf(LocalRankFilename,"neighbors.%05i.raw",rank); + NB = fopen(LocalRankFilename,"wb"); + fwrite(NeighborList,4,18*Np,NB); + fclose(NB); + + FILE *DIST; + DoubleArray DistData(7, Np); + ScaLBL_CopyToHost(DistData.data(), hq, 7*sizeof(double)*Np); + sprintf(LocalRankFilename,"h.%05i.raw",rank); + DIST = fopen(LocalRankFilename,"wb"); + fwrite(DistData.data(),8,7*Np,DIST); + fclose(DIST); + + */ FILE *OUTFILE; sprintf(LocalRankFilename,"Phase.%05i.raw",rank); @@ -940,6 +1030,17 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){ fclose(OUTFILE); DoubleArray PhaseField(Nx,Ny,Nz); + FILE *DIST; + for (int q=0; q<7; q++){ + ScaLBL_Comm->RegularLayout(Map,&hq[q*Np],PhaseField); + + sprintf(LocalRankFilename,"h%i.%05i.raw",q,rank); + DIST = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,Nx*Ny*Nz,DIST); + fclose(DIST); + + } + ScaLBL_Comm->RegularLayout(Map,Den,PhaseField); FILE *AFILE; sprintf(LocalRankFilename,"Density.%05i.raw",rank); @@ -975,7 +1076,7 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){ fwrite(PhaseField.data(),8,N,VELZ_FILE); fclose(VELZ_FILE); -/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField); + ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField); FILE *CGX_FILE; sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank); CGX_FILE = fopen(LocalRankFilename,"wb"); @@ -995,7 +1096,7 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){ CGZ_FILE = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,CGZ_FILE); fclose(CGZ_FILE); -*/ + } void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){ @@ -1031,3 +1132,151 @@ void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){ fwrite(PhaseField.data(),8,N,VELZ_FILE); fclose(VELZ_FILE); } + +void ScaLBL_FreeLeeModel::Create_DummyPhase_MGTest(){ + // Initialize communication structures in averaging domain + for (int i=0; iid[i] = Mask->id[i]; + Mask->CommInit(); + Np=Mask->PoreCount(); + //........................................................................... + if (rank==0) printf ("Create ScaLBL_Communicator \n"); + // Create a communicator for the device (will use optimized layout) + // ScaLBL_Communicator ScaLBL_Comm(Mask); // original + ScaLBL_Comm = std::shared_ptr(new ScaLBL_Communicator(Mask)); + //ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); + ScaLBL_Comm_WideHalo = std::shared_ptr(new ScaLBLWideHalo_Communicator(Mask,2)); + + // create the layout for the LBM + int Npad=(Np/16 + 2)*16; + if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); + Map.resize(Nx,Ny,Nz); Map.fill(-2); + auto neighborList= new int[18*Npad]; + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1); + comm.barrier(); + + //........................................................................... + // MAIN VARIABLES ALLOCATED HERE + //........................................................................... + // LBM variables + if (rank==0) printf ("Allocating distributions \n"); + //......................device distributions................................. + dist_mem_size = Np*sizeof(double); + neighborSize=18*(Np*sizeof(int)); + //........................................................................... + //ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); + //ScaLBL_AllocateDeviceMemory((void **) &gqbar, 19*dist_mem_size); + //ScaLBL_AllocateDeviceMemory((void **) &hq, 7*dist_mem_size); + //ScaLBL_AllocateDeviceMemory((void **) &mu_phi, dist_mem_size); + //ScaLBL_AllocateDeviceMemory((void **) &Den, dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nh); + //ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); + //ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np); + //........................................................................... + // Update GPU data structures + if (rank==0) printf ("Setting up device map and neighbor list \n"); + fflush(stdout); + int *TmpMap; + TmpMap=new int[Np]; + for (int k=1; kMap(i,j,k); + } + } + } + // check that TmpMap is valid + for (int idx=0; idxLastExterior(); idx++){ + auto n = TmpMap[idx]; + if (n > Nxh*Nyh*Nzh){ + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nxh*Nyh*Nzh-1; + } + } + for (int idx=ScaLBL_Comm->FirstInterior(); idxLastInterior(); idx++){ + auto n = TmpMap[idx]; + if ( n > Nxh*Nyh*Nzh ){ + printf("Bad value! idx=%i \n",n); + TmpMap[idx] = Nxh*Nyh*Nzh-1; + } + } + // copy the device map + ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np); + // copy the neighbor list + //ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + comm.barrier(); + + double *phase; + phase = new double[Nh]; + + for (int k=0;kid[n] + int x=i-1; + int y=j-1; + int z=k-1; + if (x<0) x=0; + if (y<0) y=0; + if (z<0) z=0; + if (x>=Nx) x=Nx-1; + if (y>=Ny) y=Ny-1; + if (z>=Nz) z=Nz-1; + int n = z*Nx*Ny+y*Nx+x; + phase[nh]=id[n]; + } + } + } + ScaLBL_CopyToDevice(Phi, phase, Nh*sizeof(double)); + ScaLBL_Comm->Barrier(); + comm.barrier(); + delete [] TmpMap; + delete [] neighborList; + delete [] phase; +} + +void ScaLBL_FreeLeeModel::MGTest(){ + + comm.barrier(); + + ScaLBL_Comm_WideHalo->Send(Phi); + ScaLBL_D3Q9_MGTest(dvcMap,Phi,ColorGrad,Nxh,Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_WideHalo->Recv(Phi); + ScaLBL_D3Q9_MGTest(dvcMap,Phi,ColorGrad,Nxh,Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np); + + //check the sum of ColorGrad + double cgx_loc = 0.0; + double cgy_loc = 0.0; + double cgz_loc = 0.0; + double cgx,cgy,cgz; + double *ColorGrad_host; + ColorGrad_host = new double [3*Np]; + ScaLBL_CopyToHost(&ColorGrad_host[0],&ColorGrad[0], 3*Np*sizeof(double)); + for (int i = ScaLBL_Comm->FirstInterior(); iLastInterior();i++){ + cgx_loc+=ColorGrad_host[0*Np+i]; + cgy_loc+=ColorGrad_host[1*Np+i]; + cgz_loc+=ColorGrad_host[2*Np+i]; + } + for (int i = 0; iLastExterior();i++){ + cgx_loc+=ColorGrad_host[0*Np+i]; + cgy_loc+=ColorGrad_host[1*Np+i]; + cgz_loc+=ColorGrad_host[2*Np+i]; + } + cgx=Dm->Comm.sumReduce( cgx_loc); + cgy=Dm->Comm.sumReduce( cgy_loc); + cgz=Dm->Comm.sumReduce( cgz_loc); + if (rank==0){ + printf("Sum of all x-component of the mixed gradient = %.2g \n",cgx); + printf("Sum of all y-component of the mixed gradient = %.2g \n",cgy); + printf("Sum of all z-component of the mixed gradient = %.2g \n",cgz); + } + + delete [] ColorGrad_host; +} diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 1e372f50..17cc6323 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -16,6 +16,9 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model #include "common/ScaLBL.h" #include "common/WideHalo.h" +#ifndef ScaLBL_FreeLeeModel_INC +#define ScaLBL_FreeLeeModel_INC + class ScaLBL_FreeLeeModel{ public: ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM); @@ -28,12 +31,17 @@ public: void ReadInput(); void Create_TwoFluid(); void Initialize_TwoFluid(); - void Run_TwoFluid(); + double Run_TwoFluid(int returntime); + void WriteDebug_TwoFluid(); void Create_SingleFluid(); void Initialize_SingleFluid(); void Run_SingleFluid(); + void WriteDebug_SingleFluid(); + // test utilities + void Create_DummyPhase_MGTest(); + void MGTest(); bool Restart,pBC; int timestep,timestepMax; @@ -73,8 +81,12 @@ public: double *Velocity; double *Pressure; + void getPhase(DoubleArray &PhaseValues); + void getPotential(DoubleArray &PressureValues, DoubleArray &MuValues); + void getVelocity(DoubleArray &Vx, DoubleArray &Vy, DoubleArray &Vz); + DoubleArray SignDist; - + private: Utilities::MPI comm; @@ -90,4 +102,4 @@ private: void AssignComponentLabels_ChemPotential_ColorGrad(); }; - +#endif diff --git a/sample_scripts/configure_ubuntu b/sample_scripts/configure_ubuntu index cccb112c..c6c3239f 100755 --- a/sample_scripts/configure_ubuntu +++ b/sample_scripts/configure_ubuntu @@ -6,7 +6,6 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \ -D CMAKE_CXX_FLAGS="-O3 -fPIC " \ -D CMAKE_CXX_STANDARD=14 \ -D MPIEXEC=mpirun \ - -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ -D CMAKE_BUILD_TYPE:STRING=Release \ -D CUDA_FLAGS="-arch sm_35" \ -D CUDA_HOST_COMPILER="/usr/bin/gcc" \ @@ -15,7 +14,7 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \ -D USE_SILO=1 \ -D SILO_LIB="/opt/arden/silo/4.10.2/lib/libsiloh5.a" \ -D SILO_DIRECTORY="/opt/arden/silo/4.10.2" \ - -D USE_NETCDF=1 \ + -D USE_NETCDF=0 \ -D NETCDF_DIRECTORY="/opt/arden/netcdf/4.6.1" \ -D USE_CUDA=0 \ -D USE_TIMER=0 \ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2c4b8431..44f869bf 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -43,6 +43,7 @@ ADD_LBPM_EXECUTABLE( TestPoissonSolver ) ADD_LBPM_EXECUTABLE( TestIonModel ) ADD_LBPM_EXECUTABLE( TestNernstPlanck ) ADD_LBPM_EXECUTABLE( TestPNP_Stokes ) +ADD_LBPM_EXECUTABLE( TestMixedGrad ) @@ -61,6 +62,7 @@ ADD_LBPM_TEST( TestMap ) ADD_LBPM_TEST( TestWideHalo ) ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db) +ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db) #ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db) ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db) ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db) diff --git a/tests/TestMixedGrad.cpp b/tests/TestMixedGrad.cpp new file mode 100644 index 00000000..6baede7c --- /dev/null +++ b/tests/TestMixedGrad.cpp @@ -0,0 +1,199 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Utilities.h" +#include "models/FreeLeeModel.h" + +inline void Initialize_Mask(ScaLBL_FreeLeeModel &LeeModel){ + // initialize a bubble + int i,j,k,n; + int rank = LeeModel.Mask->rank(); + int Nx = LeeModel.Mask->Nx; + int Ny = LeeModel.Mask->Ny; + int Nz = LeeModel.Mask->Nz; + if (rank == 0) printf(" initialize mask...\n"); + + for (k=0;kid[n]=1; + LeeModel.id[n] = LeeModel.Mask->id[n]; + } + } + } +} + + +inline void Initialize_DummyPhaseField(ScaLBL_FreeLeeModel &LeeModel, double ax, double ay, double az){ + // initialize a bubble + int i,j,k,n; + int rank = LeeModel.Mask->rank(); + int Nx = LeeModel.Mask->Nx; + int Ny = LeeModel.Mask->Ny; + int Nz = LeeModel.Mask->Nz; + if (rank == 0) printf("Setting up dummy phase field with gradient {x,y,z} = {%f , %f , %f}...\n",ax,ay,az); + + double * Dummy; + int Nh = (Nx+2)*(Ny+2)*(Nz+2); + Dummy = new double [(Nx+2)*(Ny+2)*(Nz+2)]; + for (k=0;kid[n]=1; + LeeModel.id[n] = LeeModel.Mask->id[n]; + int nh = (k+1)*(Nx+2)*(Ny+2) + (j+1)*(Nx+2) + i+1; + Dummy[nh] = ax*double(i) + ay*double(j) + az*double(k); + } + } + } + ScaLBL_CopyToDevice(LeeModel.Phi, Dummy, sizeof(double)*Nh); + + LeeModel.MGTest(); +} + +inline int MultiHaloNeighborCheck(ScaLBL_FreeLeeModel &LeeModel){ + int i,j,k,iq,stride,nread; + int Nxh = LeeModel.Nxh; + int Nyh = LeeModel.Nyh; + int Np = LeeModel.Np; + + int *TmpMap; + TmpMap = new int[Np]; + ScaLBL_CopyToHost(TmpMap, LeeModel.dvcMap, Np*sizeof(int)); + + int *neighborList; + neighborList = new int[18*Np]; + ScaLBL_CopyToHost(neighborList, LeeModel.NeighborList, 18*Np*sizeof(int)); + printf("Check stride for interior neighbors \n"); + int count = 0; + for (int n=LeeModel.ScaLBL_Comm->FirstInterior(); nLastInterior(); n++){ + // q=0 + int idx = TmpMap[n]; + k = idx/Nxh/Nyh; + j = (idx-k*Nxh*Nyh)/Nxh; + i = (idx-k*Nxh*Nyh -j*Nxh); + + // q=1 + nread = neighborList[n]; + iq = TmpMap[nread%Np]; + stride = idx - iq; + if (stride != 1){ + printf(" %i, %i, %i q = 1 stride=%i \n ",i,j,k,stride); + count++; + } + + // q=2 + nread = neighborList[n+Np]; + iq = TmpMap[nread%Np]; + stride = iq - idx; + if (stride != 1){ + printf(" %i, %i, %i q = 2 stride=%i \n ",i,j,k,stride); + count++; + } + + + // q=3 + nread = neighborList[n+2*Np]; + iq = TmpMap[nread%Np]; + stride = idx - iq; + if (stride != Nxh){ + printf(" %i, %i, %i q = 3 stride=%i \n ",i,j,k,stride); + count++; + } + + // q = 4 + nread = neighborList[n+3*Np]; + iq = TmpMap[nread%Np]; + stride = iq-idx; + if (stride != Nxh){ + printf(" %i, %i, %i q = 4 stride=%i \n ",i,j,k,stride); + count++; + } + + + // q=5 + nread = neighborList[n+4*Np]; + iq = TmpMap[nread%Np]; + stride = idx - iq; + if (stride != Nxh*Nyh){ + count++; + printf(" %i, %i, %i q = 5 stride=%i \n ",i,j,k,stride); + } + + // q = 6 + nread = neighborList[n+5*Np]; + iq = TmpMap[nread%Np]; + stride = iq - idx; + if (stride != Nxh*Nyh){ + count++; + printf(" %i, %i, %i q = 6 stride=%i \n ",i,j,k,stride); + } + + } + return count; +} + +int main( int argc, char **argv ) +{ + + // Initialize + Utilities::startup( argc, argv ); + int errors = 0; + // Load the input database + auto db = std::make_shared( argv[1] ); + + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + + if ( rank == 0 ) { + printf( "********************************************************\n" ); + printf( "Running Mixed Gradient Test \n" ); + printf( "********************************************************\n" ); + } + // Initialize compute device + int device = ScaLBL_SetDevice( rank ); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); + + PROFILE_ENABLE( 1 ); + // PROFILE_ENABLE_TRACE(); + // PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START( "Main" ); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_FreeLeeModel LeeModel( rank, nprocs, comm ); + LeeModel.ReadParams( filename ); + LeeModel.SetDomain(); + Initialize_Mask(LeeModel); + //LeeModel.Create_DummyPhase_MGTest(); + LeeModel.Create_TwoFluid(); + + errors=MultiHaloNeighborCheck(LeeModel); + + Initialize_DummyPhaseField(LeeModel,1.0, 2.0, 3.0); + LeeModel.WriteDebug_TwoFluid(); + + PROFILE_STOP( "Main" ); + PROFILE_SAVE( file, level ); + // **************************************************** + + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::shutdown(); + return errors; + +} diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 590d5b8e..d62bef0f 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -27,19 +27,24 @@ int main( int argc, char **argv ) // Initialize Utilities::startup( argc, argv ); - // Load the input database - auto db = std::make_shared( argv[1] ); - { // Limit scope so variables that contain communicators will free before MPI_Finialize Utilities::MPI comm( MPI_COMM_WORLD ); int rank = comm.getRank(); int nprocs = comm.getSize(); + std::string SimulationMode = "production"; + // Load the input database + auto db = std::make_shared( argv[1] ); + if (argc > 2) { + SimulationMode = "development"; + } if ( rank == 0 ) { printf( "********************************************************\n" ); printf( "Running Color LBM \n" ); printf( "********************************************************\n" ); + if (SimulationMode == "development") + printf("**** DEVELOPMENT MODE ENABLED *************\n"); } // Initialize compute device int device = ScaLBL_SetDevice( rank ); @@ -62,8 +67,29 @@ int main( int argc, char **argv ) ColorModel.Create(); // creating the model will create data structure to match the pore // structure and allocate variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables - ColorModel.Run(); - // ColorModel.WriteDebug(); + + if (SimulationMode == "development"){ + double MLUPS=0.0; + int timestep = 0; + int analysis_interval = ColorModel.timestepMax; + if (ColorModel.analysis_db->keyExists( "" )){ + analysis_interval = ColorModel.analysis_db->getScalar( "analysis_interval" ); + } + FlowAdaptor Adapt(ColorModel); + runAnalysis analysis(ColorModel); + while (ColorModel.timestep < ColorModel.timestepMax){ + timestep += analysis_interval; + MLUPS = ColorModel.Run(timestep); + if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); + + Adapt.MoveInterface(ColorModel); + } + } //Analysis.WriteVis(LeeModel,LeeModel.db, timestep); + + else + ColorModel.Run(); + + ColorModel.WriteDebug(); PROFILE_STOP( "Main" ); auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); diff --git a/tests/lbpm_freelee_simulator.cpp b/tests/lbpm_freelee_simulator.cpp index 3663c4e9..0f003baa 100644 --- a/tests/lbpm_freelee_simulator.cpp +++ b/tests/lbpm_freelee_simulator.cpp @@ -8,6 +8,7 @@ #include "common/Utilities.h" #include "models/FreeLeeModel.h" +#include "analysis/FreeEnergy.h" //******************************************************************* // Implementation of Free-Energy Two-Phase LBM (Lee model) @@ -52,10 +53,33 @@ int main( int argc, char **argv ) LeeModel.SetDomain(); LeeModel.ReadInput(); LeeModel.Create_TwoFluid(); + + FreeEnergyAnalyzer Analysis(LeeModel.Dm); + LeeModel.Initialize_TwoFluid(); - LeeModel.Run_TwoFluid(); - LeeModel.WriteDebug_TwoFluid(); - + + /*** RUN MAIN TIMESTEPS HERE ************/ + double MLUPS=0.0; + int timestep = 0; + int visualization_time = LeeModel.timestepMax; + if (LeeModel.vis_db->keyExists( "visualization_interval" )){ + visualization_time = LeeModel.vis_db->getScalar( "visualization_interval" ); + timestep += visualization_time; + } + while (LeeModel.timestep < LeeModel.timestepMax){ + MLUPS = LeeModel.Run_TwoFluid(timestep); + if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); + Analysis.WriteVis(LeeModel,LeeModel.db, timestep); + timestep += visualization_time; + } + //LeeModel.WriteDebug_TwoFluid(); + if (rank==0) printf("********************************************************\n"); + 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"); + // ************************************************************************ + PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_simulator" ); auto level = db->getWithDefault( "TimerLevel", 1 ); diff --git a/tests/testGlobalMassFreeLee.cpp b/tests/testGlobalMassFreeLee.cpp new file mode 100644 index 00000000..2e976854 --- /dev/null +++ b/tests/testGlobalMassFreeLee.cpp @@ -0,0 +1,101 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Utilities.h" +#include "models/FreeLeeModel.h" + +//******************************************************************* +// Implementation of Free-Energy Two-Phase LBM (Lee model) +//******************************************************************* + +int main( int argc, char **argv ) +{ + + // Initialize + Utilities::startup( argc, argv ); + + // Load the input database + auto db = std::make_shared( argv[1] ); + + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + + if (rank == 0){ + printf("********************************************************\n"); + printf("Running Free Energy Lee LBM \n"); + printf("********************************************************\n"); + } + // Initialize compute device + int device=ScaLBL_SetDevice(rank); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); + + PROFILE_ENABLE(1); + //PROFILE_ENABLE_TRACE(); + //PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START("Main"); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_FreeLeeModel LeeModel( rank,nprocs,comm ); + LeeModel.ReadParams( filename ); + LeeModel.SetDomain(); + LeeModel.ReadInput(); + LeeModel.Create_TwoFluid(); + LeeModel.Initialize_TwoFluid(); + /* check neighbors */ + + + /* Copy the initial density to test that global mass is conserved */ + int Nx = LeeModel.Dm->Nx; + int Ny = LeeModel.Dm->Ny; + int Nz = LeeModel.Dm->Nz; + DoubleArray DensityInit(Nx,Ny,Nz); + LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityInit); + + double MLUPS = LeeModel.Run_TwoFluid(LeeModel.timestepMax); + + DoubleArray DensityFinal(Nx,Ny,Nz); + LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityFinal); + + DoubleArray DensityChange(Nx,Ny,Nz); + double totalChange=0.0; + for (int k=1; k