OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'FOM' of github.com:JamesEMcClure/LBPM-WIA into FOM

Browse Source
This commit is contained in:
Mark Berrill 2021-03-26 13:21:03 -04:00
commit 651587392b
22 changed files with 3429 additions and 1666 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

181
analysis/FreeEnergy.cpp Normal file
View File

@ -0,0 +1,181 @@
#include "analysis/FreeEnergy.h"
FreeEnergyAnalyzer::FreeEnergyAnalyzer(std::shared_ptr <Domain> 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<Ion.number_ion_species; ion++){
fprintf(TIMELOG,"%.8g ",rho_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_fluctuation_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_fluctuation_global[ion]);
}
*/
fprintf(TIMELOG,"\n");
fflush(TIMELOG);
}
/* else{
fprintf(TIMELOG,"%i ",timestep);
for (int ion=0; ion<Ion.number_ion_species; ion++){
fprintf(TIMELOG,"%.8g ",rho_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_avg_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_fluctuation_local[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_fluctuation_local[ion]);
}
fflush(TIMELOG);
} */
}
void FreeEnergyAnalyzer::WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep){
auto vis_db = input_db->getDatabase( "Visualization" );
char VisName[40];
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> 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<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto VisPhase = std::make_shared<IO::Variable>();
auto VisPressure = std::make_shared<IO::Variable>();
auto VisChemicalPotential = std::make_shared<IO::Variable>();
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
if (vis_db->getWithDefault<bool>( "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<bool>( "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<bool>( "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<bool>( "save_phase", true )){
ASSERT(visData[0].vars[0]->name=="Phase");
LeeModel.getPhase(Phi);
Array<double>& PhaseData = visData[0].vars[0]->data;
fillData.copy(Phi,PhaseData);
}
if (vis_db->getWithDefault<bool>( "save_potential", true )){
ASSERT(visData[0].vars[1]->name=="Pressure");
LeeModel.getPotential(Pressure, ChemicalPotential);
Array<double>& PressureData = visData[0].vars[1]->data;
fillData.copy(Pressure,PressureData);
ASSERT(visData[0].vars[2]->name=="ChemicalPotential");
Array<double>& ChemicalPotentialData = visData[0].vars[2]->data;
fillData.copy(ChemicalPotential,ChemicalPotentialData);
}
if (vis_db->getWithDefault<bool>( "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<double>& VelxData = visData[0].vars[3]->data;
Array<double>& VelyData = visData[0].vars[4]->data;
Array<double>& 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<bool>( "write_silo", true ))
IO::writeData( timestep, visData, Dm->Comm );
/* if (vis_db->getWithDefault<bool>( "save_8bit_raw", true )){
char CurrentIDFilename[40];
sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
Averages.AggregateLabels(CurrentIDFilename);
}
*/
}

54
analysis/FreeEnergy.h Normal file
View File

@ -0,0 +1,54 @@
/*
* averaging tools for electrochemistry
*/
#ifndef FreeEnergyAnalyzer_INC
#define FreeEnergyAnalyzer_INC
#include <vector>
#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 <Domain> 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 <Domain> Dm);
~FreeEnergyAnalyzer();
void SetParams();
void Basic( ScaLBL_FreeLeeModel &LeeModel, int timestep);
void WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep);
private:
FILE *TIMELOG;
};
#endif

2
analysis/morphology.cpp
View File

@ -702,10 +702,12 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
if (rank == 0) printf(" delta=%f, growth=%f, max. displacement = %f \n",morph_delta, GrowthEstimate, MAX_DISPLACEMENT); if (rank == 0) printf(" delta=%f, growth=%f, max. displacement = %f \n",morph_delta, GrowthEstimate, MAX_DISPLACEMENT);
// Now adjust morph_delta // Now adjust morph_delta
if (fabs(GrowthEstimate - GrowthPrevious) > 0.0) {
double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious); double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious);
GrowthPrevious = GrowthEstimate; GrowthPrevious = GrowthEstimate;
morph_delta_previous = morph_delta; morph_delta_previous = morph_delta;
morph_delta += step_size; morph_delta += step_size;
}
if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0; if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0;
//MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25); //MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25);

133
analysis/runAnalysis.cpp
View File

@ -706,6 +706,139 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> 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<int>( "N_threads", 4 );
auto method = db->getWithDefault<std::string>( "load_balance", "default" );
createThreads( method, N_threads );
}
runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
/* std::shared_ptr<Database> input_db, const RankInfoStruct &rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> 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<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "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<int>( "blobid_interval" );
}
if ( db->keyExists( "visualization_interval" ) ) {
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
}
if ( db->keyExists( "subphase_analysis_interval" ) ) {
d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
}
auto restart_file = db->getScalar<std::string>( "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<IO::DomainMesh>(
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<IO::Variable>();
auto PressVar = std::make_shared<IO::Variable>();
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
auto BlobIDVar = std::make_shared<IO::Variable>();
if ( vis_db->getWithDefault<bool>( "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<bool>( "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<bool>( "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<bool>( "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<bool>( "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 // Initialize the comms
for ( int i = 0; i < 1024; i++ ) for ( int i = 0; i < 1024; i++ )
d_comm_used[i] = false; d_comm_used[i] = false;

3
analysis/runAnalysis.h
View File

@ -7,6 +7,7 @@
#include "common/Communication.h" #include "common/Communication.h"
#include "common/ScaLBL.h" #include "common/ScaLBL.h"
#include "threadpool/thread_pool.h" #include "threadpool/thread_pool.h"
#include "models/ColorModel.h"
#include <limits.h> #include <limits.h>
@ -32,6 +33,8 @@ public:
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> dm, int Np, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> dm, int Np,
bool Regular, IntArray Map ); bool Regular, IntArray Map );
runAnalysis( ScaLBL_ColorModel &ColorModel);
//! Destructor //! Destructor
~runAnalysis(); ~runAnalysis();

24
common/ScaLBL.cpp
View File

@ -1623,35 +1623,31 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component){
// Pack the distributions // Pack the distributions
//...Packing for x face(2,8,10,12,14)................................ //...Packing for x face(2,8,10,12,14)................................
ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,&Aq[Component*7*N],N);
req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag);
req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag);
//...Packing for X face(1,7,9,11,13)................................ //...Packing for X face(1,7,9,11,13)................................
ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,&Aq[Component*7*N],N);
req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag);
req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag);
//...Packing for y face(4,8,9,16,18)................................. //...Packing for y face(4,8,9,16,18).................................
ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,&Aq[Component*7*N],N);
req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y,rank_y,sendtag);
req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y,rank_Y,recvtag);
//...Packing for Y face(3,7,10,15,17)................................. //...Packing for Y face(3,7,10,15,17).................................
ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,&Aq[Component*7*N],N);
req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y,rank_Y,sendtag);
req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y,rank_y,recvtag);
//...Packing for z face(6,12,13,16,17)................................ //...Packing for z face(6,12,13,16,17)................................
ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,&Aq[Component*7*N],N);
req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z,rank_z,sendtag);
req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z,rank_Z,recvtag);
//...Packing for Z face(5,11,14,15,18)................................ //...Packing for Z face(5,11,14,15,18)................................
ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,&Aq[Component*7*N],N); ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,&Aq[Component*7*N],N);
//...................................................................................
// Send all the distributions
//...................................................................................
req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag);
req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag);
req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag);
req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag);
req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y,rank_y,sendtag);
req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y,rank_Y,recvtag);
req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y,rank_Y,sendtag);
req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y,rank_y,recvtag);
req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z,rank_z,sendtag);
req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z,rank_Z,recvtag);
req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z,rank_Z,sendtag); req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z,rank_Z,sendtag);
req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag); req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag);
} }

30
common/ScaLBL.h
View File

@ -155,6 +155,12 @@ extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *NeighborList, int *Map, double
extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi, extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi,
int start, int finish, int Np); int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, 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);
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);
extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad, int start, int finish, int Np, int Nx, int Ny, int Nz); extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *Phi, double *ColorGrad, int start, int finish, int Np, int Nx, int Ny, int Nz);
extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz); extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz);
@ -187,18 +193,25 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double
extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad,
double rhonA, double rhoB, double tauM, double W, int start, int finish, int Np); double rhonA, 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, //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); // double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, //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); // 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);
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, 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 tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np);
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); int strideY, int strideZ, int start, int finish, int Np);
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, 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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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); int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
@ -207,6 +220,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborLis
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, 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); double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np);
// BOUNDARY CONDITION ROUTINES // BOUNDARY CONDITION ROUTINES

107
common/WideHalo.cpp
View File

@ -234,59 +234,59 @@ void ScaLBLWideHalo_Communicator::Send(double *data){
//................................................................................... //...................................................................................
// Send / Recv all the phase indcator field values // Send / Recv all the phase indcator field values
//................................................................................... //...................................................................................
req1[0] = MPI_COMM_SCALBL.Isend(&sendCount_x,1,rank_x,sendtag+0); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x,sendCount_x,rank_x,sendtag+0);
req2[0] = MPI_COMM_SCALBL.Irecv(&recvCount_X,1,rank_X,recvtag+0); req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X,recvCount_X,rank_X,recvtag+0);
req1[1] = MPI_COMM_SCALBL.Isend(&sendCount_X,1,rank_X,sendtag+1); req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X,sendCount_X,rank_X,sendtag+1);
req2[1] = MPI_COMM_SCALBL.Irecv(&recvCount_x,1,rank_x,recvtag+1); req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x,recvCount_x,rank_x,recvtag+1);
req1[2] = MPI_COMM_SCALBL.Isend(&sendCount_y,1,rank_y,sendtag+2); req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y,sendCount_y,rank_y,sendtag+2);
req2[2] = MPI_COMM_SCALBL.Irecv(&recvCount_Y,1,rank_Y,recvtag+2); req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y,recvCount_Y,rank_Y,recvtag+2);
req1[3] = MPI_COMM_SCALBL.Isend(&sendCount_Y,1,rank_Y,sendtag+3); req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y,sendCount_Y,rank_Y,sendtag+3);
req2[3] = MPI_COMM_SCALBL.Irecv(&recvCount_y,1,rank_y,recvtag+3); req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y,recvCount_y,rank_y,recvtag+3);
req1[4] = MPI_COMM_SCALBL.Isend(&sendCount_z,1,rank_z,sendtag+4); req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z,sendCount_z,rank_z,sendtag+4);
req2[4] = MPI_COMM_SCALBL.Irecv(&recvCount_Z,1,rank_Z,recvtag+4); req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z,recvCount_Z,rank_Z,recvtag+4);
req1[5] = MPI_COMM_SCALBL.Isend(&sendCount_Z,1,rank_Z,sendtag+5); req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z,sendCount_Z,rank_Z,sendtag+5);
req2[5] = MPI_COMM_SCALBL.Irecv(&recvCount_z,1,rank_z,recvtag+5); req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z,recvCount_z,rank_z,recvtag+5);
req1[6] = MPI_COMM_SCALBL.Isend(&sendCount_xy,1,rank_xy,sendtag+6); req1[6] = MPI_COMM_SCALBL.Isend(sendbuf_xy,sendCount_xy,rank_xy,sendtag+6);
req2[6] = MPI_COMM_SCALBL.Irecv(&recvCount_XY,1,rank_XY,recvtag+6); req2[6] = MPI_COMM_SCALBL.Irecv(recvbuf_XY,recvCount_XY,rank_XY,recvtag+6);
req1[7] = MPI_COMM_SCALBL.Isend(&sendCount_XY,1,rank_XY,sendtag+7); req1[7] = MPI_COMM_SCALBL.Isend(sendbuf_XY,sendCount_XY,rank_XY,sendtag+7);
req2[7] = MPI_COMM_SCALBL.Irecv(&recvCount_xy,1,rank_xy,recvtag+7); req2[7] = MPI_COMM_SCALBL.Irecv(recvbuf_xy,recvCount_xy,rank_xy,recvtag+7);
req1[8] = MPI_COMM_SCALBL.Isend(&sendCount_Xy,1,rank_Xy,sendtag+8); req1[8] = MPI_COMM_SCALBL.Isend(sendbuf_Xy,sendCount_Xy,rank_Xy,sendtag+8);
req2[8] = MPI_COMM_SCALBL.Irecv(&recvCount_xY,1,rank_xY,recvtag+8); req2[8] = MPI_COMM_SCALBL.Irecv(recvbuf_xY,recvCount_xY,rank_xY,recvtag+8);
req1[9] = MPI_COMM_SCALBL.Isend(&sendCount_xY,1,rank_xY,sendtag+9); req1[9] = MPI_COMM_SCALBL.Isend(sendbuf_xY,sendCount_xY,rank_xY,sendtag+9);
req2[9] = MPI_COMM_SCALBL.Irecv(&recvCount_Xy,1,rank_Xy,recvtag+9); req2[9] = MPI_COMM_SCALBL.Irecv(recvbuf_Xy,recvCount_Xy,rank_Xy,recvtag+9);
req1[10] = MPI_COMM_SCALBL.Isend(&sendCount_xz,1,rank_xz,sendtag+10); req1[10] = MPI_COMM_SCALBL.Isend(sendbuf_xz,sendCount_xz,rank_xz,sendtag+10);
req2[10] = MPI_COMM_SCALBL.Irecv(&recvCount_XZ,1,rank_XZ,recvtag+10); req2[10] = MPI_COMM_SCALBL.Irecv(recvbuf_XZ,recvCount_XZ,rank_XZ,recvtag+10);
req1[11] = MPI_COMM_SCALBL.Isend(&sendCount_XZ,1,rank_XZ,sendtag+11); req1[11] = MPI_COMM_SCALBL.Isend(sendbuf_XZ,sendCount_XZ,rank_XZ,sendtag+11);
req2[11] = MPI_COMM_SCALBL.Irecv(&recvCount_xz,1,rank_xz,recvtag+11); req2[11] = MPI_COMM_SCALBL.Irecv(recvbuf_xz,recvCount_xz,rank_xz,recvtag+11);
req1[12] = MPI_COMM_SCALBL.Isend(&sendCount_Xz,1,rank_Xz,sendtag+12); req1[12] = MPI_COMM_SCALBL.Isend(sendbuf_Xz,sendCount_Xz,rank_Xz,sendtag+12);
req2[12] = MPI_COMM_SCALBL.Irecv(&recvCount_xZ,1,rank_xZ,recvtag+12); req2[12] = MPI_COMM_SCALBL.Irecv(recvbuf_xZ,recvCount_xZ,rank_xZ,recvtag+12);
req1[13] = MPI_COMM_SCALBL.Isend(&sendCount_xZ,1,rank_xZ,sendtag+13); req1[13] = MPI_COMM_SCALBL.Isend(sendbuf_xZ,sendCount_xZ,rank_xZ,sendtag+13);
req2[13] = MPI_COMM_SCALBL.Irecv(&recvCount_Xz,1,rank_Xz,recvtag+13); req2[13] = MPI_COMM_SCALBL.Irecv(recvbuf_Xz,recvCount_Xz,rank_Xz,recvtag+13);
req1[14] = MPI_COMM_SCALBL.Isend(&sendCount_yz,1,rank_yz,sendtag+14); req1[14] = MPI_COMM_SCALBL.Isend(sendbuf_yz,sendCount_yz,rank_yz,sendtag+14);
req2[14] = MPI_COMM_SCALBL.Irecv(&recvCount_YZ,1,rank_YZ,recvtag+14); req2[14] = MPI_COMM_SCALBL.Irecv(recvbuf_YZ,recvCount_YZ,rank_YZ,recvtag+14);
req1[15] = MPI_COMM_SCALBL.Isend(&sendCount_YZ,1,rank_YZ,sendtag+15); req1[15] = MPI_COMM_SCALBL.Isend(sendbuf_YZ,sendCount_YZ,rank_YZ,sendtag+15);
req2[15] = MPI_COMM_SCALBL.Irecv(&recvCount_yz,1,rank_yz,recvtag+15); req2[15] = MPI_COMM_SCALBL.Irecv(recvbuf_yz,recvCount_yz,rank_yz,recvtag+15);
req1[16] = MPI_COMM_SCALBL.Isend(&sendCount_Yz,1,rank_Yz,sendtag+16); req1[16] = MPI_COMM_SCALBL.Isend(sendbuf_Yz,sendCount_Yz,rank_Yz,sendtag+16);
req2[16] = MPI_COMM_SCALBL.Irecv(&recvCount_yZ,1,rank_yZ,recvtag+16); req2[16] = MPI_COMM_SCALBL.Irecv(recvbuf_yZ,recvCount_yZ,rank_yZ,recvtag+16);
req1[17] = MPI_COMM_SCALBL.Isend(&sendCount_yZ,1,rank_yZ,sendtag+17); req1[17] = MPI_COMM_SCALBL.Isend(sendbuf_yZ,sendCount_yZ,rank_yZ,sendtag+17);
req2[17] = MPI_COMM_SCALBL.Irecv(&recvCount_Yz,1,rank_Yz,recvtag+17); req2[17] = MPI_COMM_SCALBL.Irecv(recvbuf_Yz,recvCount_Yz,rank_Yz,recvtag+17);
/* Corners */ /* Corners */
req1[18] = MPI_COMM_SCALBL.Isend(&sendCount_xyz,1,rank_xyz,sendtag+18); req1[18] = MPI_COMM_SCALBL.Isend(sendbuf_xyz,sendCount_xyz,rank_xyz,sendtag+18);
req2[18] = MPI_COMM_SCALBL.Irecv(&recvCount_XYZ,1,rank_XYZ,recvtag+18); req2[18] = MPI_COMM_SCALBL.Irecv(recvbuf_XYZ,recvCount_XYZ,rank_XYZ,recvtag+18);
req1[19] = MPI_COMM_SCALBL.Isend(&sendCount_XYz,1,rank_XYz,sendtag+19); req1[19] = MPI_COMM_SCALBL.Isend(sendbuf_XYz,sendCount_XYz,rank_XYz,sendtag+19);
req2[19] = MPI_COMM_SCALBL.Irecv(&recvCount_xyZ,1,rank_xyZ,recvtag+19); req2[19] = MPI_COMM_SCALBL.Irecv(recvbuf_xyZ,recvCount_xyZ,rank_xyZ,recvtag+19);
req1[20] = MPI_COMM_SCALBL.Isend(&sendCount_Xyz,1,rank_Xyz,sendtag+20); req1[20] = MPI_COMM_SCALBL.Isend(sendbuf_Xyz,sendCount_Xyz,rank_Xyz,sendtag+20);
req2[20] = MPI_COMM_SCALBL.Irecv(&recvCount_xYZ,1,rank_xYZ,recvtag+20); req2[20] = MPI_COMM_SCALBL.Irecv(recvbuf_xYZ,recvCount_xYZ,rank_xYZ,recvtag+20);
req1[21] = MPI_COMM_SCALBL.Isend(&sendCount_xYz,1,rank_xYz,sendtag+21); req1[21] = MPI_COMM_SCALBL.Isend(sendbuf_xYz,sendCount_xYz,rank_xYz,sendtag+21);
req2[21] = MPI_COMM_SCALBL.Irecv(&recvCount_XyZ,1,rank_XyZ,recvtag+21); req2[21] = MPI_COMM_SCALBL.Irecv(recvbuf_XyZ,recvCount_XyZ,rank_XyZ,recvtag+21);
req1[22] = MPI_COMM_SCALBL.Isend(&sendCount_xyZ,1,rank_xyZ,sendtag+22); req1[22] = MPI_COMM_SCALBL.Isend(sendbuf_xyZ,sendCount_xyZ,rank_xyZ,sendtag+22);
req2[22] = MPI_COMM_SCALBL.Irecv(&recvCount_XYz,1,rank_XYz,recvtag+22); req2[22] = MPI_COMM_SCALBL.Irecv(recvbuf_XYz,recvCount_XYz,rank_XYz,recvtag+22);
req1[23] = MPI_COMM_SCALBL.Isend(&sendCount_XYZ,1,rank_XYZ,sendtag+23); req1[23] = MPI_COMM_SCALBL.Isend(sendbuf_XYZ,sendCount_XYZ,rank_XYZ,sendtag+23);
req2[23] = MPI_COMM_SCALBL.Irecv(&recvCount_xyz,1,rank_xyz,recvtag+23); req2[23] = MPI_COMM_SCALBL.Irecv(recvbuf_xyz,recvCount_xyz,rank_xyz,recvtag+23);
req1[24] = MPI_COMM_SCALBL.Isend(&sendCount_XyZ,1,rank_XyZ,sendtag+24); req1[24] = MPI_COMM_SCALBL.Isend(sendbuf_XyZ,sendCount_XyZ,rank_XyZ,sendtag+24);
req2[24] = MPI_COMM_SCALBL.Irecv(&recvCount_xYz,1,rank_xYz,recvtag+24); req2[24] = MPI_COMM_SCALBL.Irecv(recvbuf_xYz,recvCount_xYz,rank_xYz,recvtag+24);
req1[25] = MPI_COMM_SCALBL.Isend(&sendCount_xYZ,1,rank_xYZ,sendtag+25); req1[25] = MPI_COMM_SCALBL.Isend(sendbuf_xYZ,sendCount_xYZ,rank_xYZ,sendtag+25);
req2[25] = MPI_COMM_SCALBL.Irecv(&recvCount_Xyz,1,rank_Xyz,recvtag+25); req2[25] = MPI_COMM_SCALBL.Irecv(recvbuf_Xyz,recvCount_Xyz,rank_Xyz,recvtag+25);
//................................................................................... //...................................................................................
} }
@ -302,6 +302,9 @@ void ScaLBLWideHalo_Communicator::Recv(double *data){
Utilities::MPI::waitAll(26,req2); Utilities::MPI::waitAll(26,req2);
ScaLBL_DeviceBarrier(); ScaLBL_DeviceBarrier();
//................................................................................... //...................................................................................
//printf("Ready to unpack %i to x\n",recvCount_x);
//printf(" print first 10 values...\n");
//for (int idx=0; idx<10; idx++) printf(" recvBuf[%i]=%f \n",idx,recvbuf_x[idx]);
ScaLBL_Scalar_Unpack(dvcRecvList_x, recvCount_x,recvbuf_x, data, Nh); ScaLBL_Scalar_Unpack(dvcRecvList_x, recvCount_x,recvbuf_x, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_y, recvCount_y,recvbuf_y, data, Nh); ScaLBL_Scalar_Unpack(dvcRecvList_y, recvCount_y,recvbuf_y, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_X, recvCount_X,recvbuf_X, data, Nh); ScaLBL_Scalar_Unpack(dvcRecvList_X, recvCount_X,recvbuf_X, data, Nh);

194
cpu/Color.cpp
View File

@ -2489,10 +2489,200 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
} }
} }
extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, 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){
int nr1,nr2,nr3,nr4,nr5,nr6;
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; n<finish; n++){
/* neighbors */
nr1 = neighborList[n+0*Np];
nr2 = neighborList[n+1*Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
/* load velocity */
ux = Vel[n];
uy = Vel[Np+n];
uz = Vel[2*Np+n];
/* load color gradient */
nx = ColorGrad[n];
ny = ColorGrad[Np+n];
nz = ColorGrad[2*Np+n];
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) 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; n<finish; n++){
/* load velocity */
ux = Vel[n];
uy = Vel[Np+n];
uz = Vel[2*Np+n];
/* load color gradient */
nx = ColorGrad[n];
ny = ColorGrad[Np+n];
nz = ColorGrad[2*Np+n];
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) 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, 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){ double *Den, double *Phi, int start, int finish, int Np){
int idx,n,nread; int idx,nread;
double fq,nA,nB; double fq,nA,nB;
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){
@ -2579,7 +2769,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double
extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi, extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq, double *Den, double *Phi,
int start, int finish, int Np){ int start, int finish, int Np){
int idx,n,nread; int idx,nread;
double fq,nA,nB; double fq,nA,nB;
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){

520
cpu/FreeLee.cpp
View File

@ -1,4 +1,5 @@
#include <math.h> #include <math.h>
#include <stdio.h>
#define STOKES #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[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; 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, extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad,
@ -101,6 +104,7 @@ extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, doubl
nz = nz/ColorMag_temp; nz = nz/ColorMag_temp;
theta = M*cs2_inv*(1-4.0*phi*phi)/W; 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[0*Np+idx]=0.3333333333333333*(phi);
hq[1*Np+idx]=0.1111111111111111*(phi+theta*nx); hq[1*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, 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){ double rhoA, double rhoB, int start, int finish, int Np){
int idx,n,nread; int idx,nread;
double fq,phi; double fq,phi;
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){
@ -161,12 +165,15 @@ extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int
// save the phase indicator field // save the phase indicator field
idx = Map[n]; idx = Map[n];
Phi[idx] = phi; Phi[idx] = phi;
} }
} }
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, 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){ double rhoA, double rhoB, int start, int finish, int Np){
int idx,n; int idx;
double fq,phi; double fq,phi;
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){
@ -207,11 +214,189 @@ extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq,
} }
} }
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, 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 tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int 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<finish; n++){
/* load phase indicator field */
idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=1
nr1 = neighborList[n];
nr2 = neighborList[n+Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
//q=0
h0 = hq[n];
//q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
//-------------------------------- BGK collison for phase field ---------------------------------//
h0 -= (h0 - 0.3333333333333333*phi)/tauM;
h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[nr2] = h1;
hq[nr1] = h2;
hq[nr4] = h3;
hq[nr3] = h4;
hq[nr6] = h5;
hq[nr5] = h6;
//........................................................................
}
}
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){
int idx,n;
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<finish; n++){
/* load phase indicator field */
idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
h0 = hq[n];
h1 = hq[2*Np+n];
h2 = hq[Np+n];
h3 = hq[4*Np+n];
h4 = hq[3*Np+n];
h5 = hq[6*Np+n];
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
h0 -= (h0 - 0.3333333333333333*phi)/tauM;
h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[Np+n] = h1;
hq[2*Np+n] = h2;
hq[3*Np+n] = h3;
hq[4*Np+n] = h4;
hq[5*Np+n] = h5;
hq[6*Np+n] = h6;
//........................................................................
}
}
extern "C" void ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
int idx,n;
double h0,h1,h2,h3,h4,h5,h6;
double phi;
for (int n=start; n<finish; n++){
h0 = hq[n];
h1 = hq[1*Np+n];
h2 = hq[2*Np+n];
h3 = hq[3*Np+n];
h4 = hq[4*Np+n];
h5 = hq[5*Np+n];
h6 = hq[6*Np+n];
phi = h0+h1+h2+h3+h4+h5+h6;
// save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// save the phase indicator field
idx = Map[n];
Phi[idx] = phi;
}
}
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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int nn,nn2x,ijk;
int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18; int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
@ -228,20 +413,22 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, dou
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 // double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -383,9 +570,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, dou
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -741,62 +928,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, dou
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -813,11 +944,11 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, dou
} }
} }
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, 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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int nn,nn2x,ijk;
//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18; //int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
@ -834,20 +965,22 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 //double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
for (int n=start; n<finish; n++){ for (int n=start; n<finish; n++){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -989,9 +1122,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -1053,6 +1186,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0); ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0); uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0); uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure //compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17) p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz); +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
@ -1330,62 +1464,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[2*Np+n];
//q=2
h2 = hq[1*Np+n];
//q=3
h3 = hq[4*Np+n];
//q=4
h4 = hq[3*Np+n];
//q=5
h5 = hq[6*Np+n];
//q=6
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -1405,7 +1483,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, 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){ double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
int n;
int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18; int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
@ -1672,7 +1749,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborLis
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, 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){ double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){
int n;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
// distribution functions // distribution functions
@ -1916,3 +1992,167 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, d
Pressure[n] = p; Pressure[n] = p;
} }
} }
extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad,int strideY, int strideZ, int start, int finish, int Np){
int nn,nn2x,ijk;
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m3,m5,m7;
double mm1,mm2,mm4,mm6,mm8,mm9,mm10,mm11,mm12,mm13,mm14,mm15,mm16,mm17,mm18;
double mm3,mm5,mm7;
//double nx,ny,nz;//normal color gradient
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
double phi;
for (int n=start; n<finish; n++){
// Get the 1D index based on regular data layout
ijk = Map[n];
phi = Phi[ijk];// load phase field
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
// compute mixed difference (Eq.30, A.Fukhari et al. JCP 315(2016) 434-457)
//........................................................................
nn2x = ijk-2; // neighbor index (get convention)
mm1 = Phi[nn2x]; // get neighbor for phi - 1
mm1 = 0.25*(-mm1+5.0*m1-3.0*phi-m2);
//........................................................................
nn2x = ijk+2; // neighbor index (get convention)
mm2 = Phi[nn2x]; // get neighbor for phi - 2
mm2 = 0.25*(-mm2+5.0*m2-3.0*phi-m1);
//........................................................................
nn2x = ijk-strideY*2; // neighbor index (get convention)
mm3 = Phi[nn2x]; // get neighbor for phi - 3
mm3 = 0.25*(-mm3+5.0*m3-3.0*phi-m4);
//........................................................................
nn2x = ijk+strideY*2; // neighbor index (get convention)
mm4 = Phi[nn2x]; // get neighbor for phi - 4
mm4 = 0.25*(-mm4+5.0*m4-3.0*phi-m3);
//........................................................................
nn2x = ijk-strideZ*2; // neighbor index (get convention)
mm5 = Phi[nn2x]; // get neighbor for phi - 5
mm5 = 0.25*(-mm5+5.0*m5-3.0*phi-m6);
//........................................................................
nn2x = ijk+strideZ*2; // neighbor index (get convention)
mm6 = Phi[nn2x]; // get neighbor for phi - 6
mm6 = 0.25*(-mm6+5.0*m6-3.0*phi-m5);
//........................................................................
nn2x = ijk-strideY*2-2; // neighbor index (get convention)
mm7 = Phi[nn2x]; // get neighbor for phi - 7
mm7 = 0.25*(-mm7+5.0*m7-3.0*phi-m8);
//........................................................................
nn2x = ijk+strideY*2+2; // neighbor index (get convention)
mm8 = Phi[nn2x]; // get neighbor for phi - 8
mm8 = 0.25*(-mm8+5.0*m8-3.0*phi-m7);
//........................................................................
nn2x = ijk+strideY*2-2; // neighbor index (get convention)
mm9 = Phi[nn2x]; // get neighbor for phi - 9
mm9 = 0.25*(-mm9+5.0*m9-3.0*phi-m10);
//........................................................................
nn2x = ijk-strideY*2+2; // neighbor index (get convention)
mm10 = Phi[nn2x]; // get neighbor for phi - 10
mm10 = 0.25*(-mm10+5.0*m10-3.0*phi-m9);
//........................................................................
nn2x = ijk-strideZ*2-2; // neighbor index (get convention)
mm11 = Phi[nn2x]; // get neighbor for phi - 11
mm11 = 0.25*(-mm11+5.0*m11-3.0*phi-m12);
//........................................................................
nn2x = ijk+strideZ*2+2; // neighbor index (get convention)
mm12 = Phi[nn2x]; // get neighbor for phi - 12
mm12 = 0.25*(-mm12+5.0*m12-3.0*phi-m11);
//........................................................................
nn2x = ijk+strideZ*2-2; // neighbor index (get convention)
mm13 = Phi[nn2x]; // get neighbor for phi - 13
mm13 = 0.25*(-mm13+5.0*m13-3.0*phi-m14);
//........................................................................
nn2x = ijk-strideZ*2+2; // neighbor index (get convention)
mm14 = Phi[nn2x]; // get neighbor for phi - 14
mm14 = 0.25*(-mm14+5.0*m14-3.0*phi-m13);
//........................................................................
nn2x = ijk-strideZ*2-strideY*2; // neighbor index (get convention)
mm15 = Phi[nn2x]; // get neighbor for phi - 15
mm15 = 0.25*(-mm15+5.0*m15-3.0*phi-m16);
//........................................................................
nn2x = ijk+strideZ*2+strideY*2; // neighbor index (get convention)
mm16 = Phi[nn2x]; // get neighbor for phi - 16
mm16 = 0.25*(-mm16+5.0*m16-3.0*phi-m15);
//........................................................................
nn2x = ijk+strideZ*2-strideY*2; // neighbor index (get convention)
mm17 = Phi[nn2x]; // get neighbor for phi - 17
mm17 = 0.25*(-mm17+5.0*m17-3.0*phi-m18);
//........................................................................
nn2x = ijk-strideZ*2+strideY*2; // neighbor index (get convention)
mm18 = Phi[nn2x]; // get neighbor for phi - 18
mm18 = 0.25*(-mm18+5.0*m18-3.0*phi-m17);
//............Compute the Color Gradient...................................
//nx = -3.0*1.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
//ny = -3.0*1.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
//nz = -3.0*1.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
ColorGrad[0*Np+n] = mgx;
ColorGrad[1*Np+n] = mgy;
ColorGrad[2*Np+n] = mgz;
}
}

454
cuda/FreeLee.cu
View File

@ -1,4 +1,6 @@
#include <math.h> #include <math.h>
#include <stdio.h>
#include <cuda_profiler_api.h>
#define STOKES #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, __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, int start, int finish, int Np){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
double fq,phi; 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<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -196,33 +204,168 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start; n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){ if ( n<finish ){
// q=0
fq = hq[n]; /* load phase indicator field */
phi = fq; idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=1 // q=1
fq = hq[2*Np+n]; nr1 = neighborList[n];
phi += fq; nr2 = neighborList[n+Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
// f2 = hq[10*Np+n]; //q=0
fq = hq[1*Np+n]; h0 = hq[n];
phi += fq; //q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
// q=3 //-------------------------------- BGK collison for phase field ---------------------------------//
fq = hq[4*Np+n]; h0 -= (h0 - 0.3333333333333333*phi)/tauM;
phi += fq; h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
// q = 4 /*Update the distributions */
fq = hq[3*Np+n]; // q = 0
phi += fq; hq[n] = h0;
hq[nr2] = h1;
hq[nr1] = h2;
hq[nr4] = h3;
hq[nr3] = h4;
hq[nr6] = h5;
hq[nr5] = h6;
//........................................................................
// q=5 }
fq = hq[6*Np+n]; }
phi += fq; }
// q = 6
fq = hq[5*Np+n]; __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
phi += fq; double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
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;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
/* load phase indicator field */
idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
h0 = hq[n];
h1 = hq[2*Np+n];
h2 = hq[Np+n];
h3 = hq[4*Np+n];
h4 = hq[3*Np+n];
h5 = hq[6*Np+n];
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
h0 -= (h0 - 0.3333333333333333*phi)/tauM;
h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[Np+n] = h1;
hq[2*Np+n] = h2;
hq[3*Np+n] = h3;
hq[4*Np+n] = h4;
hq[5*Np+n] = h5;
hq[6*Np+n] = h6;
//........................................................................
}
}
}
__global__ void dvc_ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
int idx,n;
double h0,h1,h2,h3,h4,h5,h6;
double phi;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
h0 = hq[n];
h1 = hq[1*Np+n];
h2 = hq[2*Np+n];
h3 = hq[3*Np+n];
h4 = hq[4*Np+n];
h5 = hq[5*Np+n];
h6 = hq[6*Np+n];
phi = h0+h1+h2+h3+h4+h5+h6;
// save the number densities // save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
@ -234,8 +377,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
} }
} }
__global__ void dvc_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, __global__ void dvc_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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
@ -255,11 +398,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 // double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -268,12 +410,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -415,9 +559,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -773,62 +917,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -845,12 +933,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
} }
} }
__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, __global__ void dvc_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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
double chem;//chemical potential double chem;//chemical potential
@ -866,10 +953,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 //double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
@ -879,12 +966,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -1026,9 +1115,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -1090,6 +1179,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0); ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0); uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0); uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure //compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17) p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz); +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
@ -1367,62 +1457,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[2*Np+n];
//q=2
h2 = hq[1*Np+n];
//q=3
h3 = hq[4*Np+n];
//q=4
h4 = hq[3*Np+n];
//q=5
h5 = hq[6*Np+n];
//q=6
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -1435,9 +1469,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ColorGrad[0*Np+n] = nx; ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny; ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz; ColorGrad[2*Np+n] = nz;
}
} }
}
} }
__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
@ -1969,48 +2003,120 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dis
extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init<<<NBLOCKS,NTHREADS >>>( 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){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<<NBLOCKS,NTHREADS >>>( 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, 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){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
} dvc_ScaLBL_FreeLeeModel_PhaseField_Init<<<NBLOCKS,NTHREADS >>>(Map, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, start, finish, Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, cudaError_t err = cudaGetLastError();
double rhoA, double rhoB, int start, int finish, int Np){ if (cudaSuccess != err){
printf("CUDA error in ScaLBL_FreeLeeModel_PhaseField_Init: %s \n",cudaGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np)
double rhoA, double rhoB, int start, int finish, int Np){ {
cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField, cudaFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField<<<NBLOCKS,NTHREADS >>>(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, 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 tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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<<<NBLOCKS,NTHREADS >>>( 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<<<NBLOCKS,NTHREADS >>>( 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){ int strideY, int strideZ, int start, int finish, int Np){
cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel, cudaFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(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 *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
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 kappa, double beta, double W, double Fx, double Fy, double Fz,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel, cudaFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(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, 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){ 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<<<NBLOCKS,NTHREADS >>>(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, 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){ 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<<<NBLOCKS,NTHREADS >>>(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){
}

457
hip/FreeLee.cu
View File

@ -1,11 +1,12 @@
#include <math.h> #include <math.h>
#include <stdio.h>
#include "hip/hip_runtime.h" #include "hip/hip_runtime.h"
#define STOKES
#define NBLOCKS 1024 #define NBLOCKS 1024
#define NTHREADS 256 #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) __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; 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, __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, int start, int finish, int Np){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
double fq,phi; 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<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -197,33 +204,168 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start; n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){ if ( n<finish ){
// q=0
fq = hq[n]; /* load phase indicator field */
phi = fq; idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=1 // q=1
fq = hq[2*Np+n]; nr1 = neighborList[n];
phi += fq; nr2 = neighborList[n+Np];
nr3 = neighborList[n+2*Np];
nr4 = neighborList[n+3*Np];
nr5 = neighborList[n+4*Np];
nr6 = neighborList[n+5*Np];
// f2 = hq[10*Np+n]; //q=0
fq = hq[1*Np+n]; h0 = hq[n];
phi += fq; //q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
// q=3 //-------------------------------- BGK collison for phase field ---------------------------------//
fq = hq[4*Np+n]; h0 -= (h0 - 0.3333333333333333*phi)/tauM;
phi += fq; h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
// q = 4 /*Update the distributions */
fq = hq[3*Np+n]; // q = 0
phi += fq; hq[n] = h0;
hq[nr2] = h1;
hq[nr1] = h2;
hq[nr4] = h3;
hq[nr3] = h4;
hq[nr6] = h5;
hq[nr5] = h6;
//........................................................................
// q=5 }
fq = hq[6*Np+n]; }
phi += fq; }
// q = 6
fq = hq[5*Np+n]; __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLee_PhaseField( int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
phi += fq; double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
int idx,n;
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;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
/* load phase indicator field */
idx = Map[n];
phi = Phi[idx];
/* velocity */
ux = Vel[0*Np+n];
uy = Vel[1*Np+n];
uz = Vel[2*Np+n];
/*color gradient */
nx = ColorGrad[0*Np+n];
ny = ColorGrad[1*Np+n];
nz = ColorGrad[2*Np+n];
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C < 1.0e-12) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
h0 = hq[n];
h1 = hq[2*Np+n];
h2 = hq[Np+n];
h3 = hq[4*Np+n];
h4 = hq[3*Np+n];
h5 = hq[6*Np+n];
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
h0 -= (h0 - 0.3333333333333333*phi)/tauM;
h1 -= (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h2 -= (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - factor*phi*phi))/W)/tauM;
h3 -= (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h4 -= (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - factor*phi*phi))/W)/tauM;
h5 -= (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
h6 -= (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - factor*phi*phi))/W)/tauM;
//........................................................................
/*Update the distributions */
// q = 0
hq[n] = h0;
hq[Np+n] = h1;
hq[2*Np+n] = h2;
hq[3*Np+n] = h3;
hq[4*Np+n] = h4;
hq[5*Np+n] = h5;
hq[6*Np+n] = h6;
//........................................................................
}
}
}
__global__ void dvc_ScaLBL_D3Q7_ComputePhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, int start, int finish, int Np){
int idx,n;
double h0,h1,h2,h3,h4,h5,h6;
double phi;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if ( n<finish ){
h0 = hq[n];
h1 = hq[1*Np+n];
h2 = hq[2*Np+n];
h3 = hq[3*Np+n];
h4 = hq[4*Np+n];
h5 = hq[5*Np+n];
h6 = hq[6*Np+n];
phi = h0+h1+h2+h3+h4+h5+h6;
// save the number densities // save the number densities
Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); Den[n] = rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
@ -235,8 +377,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double
} }
} }
__global__ void dvc_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, __global__ void dvc_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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
@ -256,11 +398,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 // double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){ for (int s=0; s<S; s++){
@ -269,12 +410,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -416,9 +559,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -774,62 +917,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[nr1];
//q=2
h2 = hq[nr2];
//q=3
h3 = hq[nr3];
//q=4
h4 = hq[nr4];
//q=5
h5 = hq[nr5];
//q=6
h6 = hq[nr6];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -846,12 +933,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map,
} }
} }
__global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad, __global__ void dvc_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 tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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){ int strideY, int strideZ, int start, int finish, int Np){
int n,nn,nn2x,ijk; int n,nn,nn2x,ijk;
//int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
double ux,uy,uz;//fluid velocity double ux,uy,uz;//fluid velocity
double p;//pressure double p;//pressure
double chem;//chemical potential double chem;//chemical potential
@ -867,10 +953,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor double mgx,mgy,mgz;//mixed gradient reaching secondary neighbor
//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; //double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field //double h0,h1,h2,h3,h4,h5,h6;//distributions for LB phase field
double tau;//position dependent LB relaxation time for fluid double tau;//position dependent LB relaxation time for fluid
double C,theta; //double C,theta;
double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7 //double M = 2.0/9.0*(tauM-0.5);//diffusivity (or mobility) for the phase field D3Q7
// for (int n=start; n<finish; n++){ // for (int n=start; n<finish; n++){
int S = Np/NBLOCKS/NTHREADS + 1; int S = Np/NBLOCKS/NTHREADS + 1;
@ -880,12 +966,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
if ( n<finish ){ if ( n<finish ){
rho0 = Den[n];//load density rho0 = Den[n];//load density
phi = Phi[n];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// Get the 1D index based on regular data layout // Get the 1D index based on regular data layout
ijk = Map[n]; ijk = Map[n];
phi = Phi[ijk];// load phase field
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
// COMPUTE THE COLOR GRADIENT // COMPUTE THE COLOR GRADIENT
//........................................................................ //........................................................................
//.................Read Phase Indicator Values............................ //.................Read Phase Indicator Values............................
@ -1027,9 +1115,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian chem = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*phi));//intermediate var, i.e. the laplacian
chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem; chem = 4.0*beta*phi*(phi+1.0)*(phi-1.0)-kappa*chem;
//............Compute the Mixed Gradient................................... //............Compute the Mixed Gradient...................................
mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14))*0.25;//the factor of 0.25 comes from the denominator of Eq.30 mgx = -3.0*1.0/18.0*(mm1-mm2+0.5*(mm7-mm8+mm9-mm10+mm11-mm12+mm13-mm14));
mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18))*0.25; mgy = -3.0*1.0/18.0*(mm3-mm4+0.5*(mm7-mm8-mm9+mm10+mm15-mm16+mm17-mm18));
mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18))*0.25; mgz = -3.0*1.0/18.0*(mm5-mm6+0.5*(mm11-mm12-mm13+mm14+mm15-mm16-mm17+mm18));
// q=0 // q=0
m0 = dist[n]; m0 = dist[n];
@ -1091,6 +1179,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0); ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(chem*nx+Fx)/3.0);
uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0); uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(chem*ny+Fy)/3.0);
uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0); uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(chem*nz+Fz)/3.0);
//compute pressure //compute pressure
p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17) p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17)
+0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz); +0.5*(rhoA-rhoB)/2.0/3.0*(ux*nx+uy*ny+uz*nz);
@ -1368,62 +1457,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)))); 0.1111111111111111*(-4*chem + (rhoA - rhoB)*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))));
//----------------------------------------------------------------------------------------------------------------------------------------// //----------------------------------------------------------------------------------------------------------------------------------------//
// ----------------------------- compute phase field evolution ----------------------------------------
//Normalize the Color Gradient
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//compute surface tension-related parameter
theta = M*4.5*(1-4.0*phi*phi)/W;
//load distributions of phase field
//q=0
h0 = hq[n];
//q=1
h1 = hq[2*Np+n];
//q=2
h2 = hq[1*Np+n];
//q=3
h3 = hq[4*Np+n];
//q=4
h4 = hq[3*Np+n];
//q=5
h5 = hq[6*Np+n];
//q=6
h6 = hq[5*Np+n];
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
// q = 1
hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 2
hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
// q = 3
hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 4
hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
// q = 5
hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
// q = 6
hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
//........................................................................
//Update velocity on device //Update velocity on device
Vel[0*Np+n] = ux; Vel[0*Np+n] = ux;
Vel[1*Np+n] = uy; Vel[1*Np+n] = uy;
@ -1436,9 +1469,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, dou
ColorGrad[0*Np+n] = nx; ColorGrad[0*Np+n] = nx;
ColorGrad[1*Np+n] = ny; ColorGrad[1*Np+n] = ny;
ColorGrad[2*Np+n] = nz; ColorGrad[2*Np+n] = nz;
}
} }
}
} }
__global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, __global__ void dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
@ -1970,48 +2003,120 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dis
extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init<<<NBLOCKS,NTHREADS >>>( 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){ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np){
dvc_ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init<<<NBLOCKS,NTHREADS >>>( 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, 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){ double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np){
} dvc_ScaLBL_FreeLeeModel_PhaseField_Init<<<NBLOCKS,NTHREADS >>>(Map, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, start, finish, Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, hipError_t err = hipGetLastError();
double rhoA, double rhoB, int start, int finish, int Np){ if (hipSuccess != err){
printf("CUDA error in ScaLBL_FreeLeeModel_PhaseField_Init: %s \n",hipGetErrorString(err));
}
} }
extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi, double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np)
double rhoA, double rhoB, int start, int finish, int Np){ {
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q7_AAodd_FreeLee_PhaseField<<<NBLOCKS,NTHREADS >>>(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, 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 tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, 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<<<NBLOCKS,NTHREADS >>>( 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<<<NBLOCKS,NTHREADS >>>( 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){ int strideY, int strideZ, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAodd_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(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 *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
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 kappa, double beta, double W, double Fx, double Fy, double Fz,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np){ int strideY, int strideZ, int start, int finish, int Np){
hipFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel, hipFuncCachePreferL1);
dvc_ScaLBL_D3Q19_AAeven_FreeLeeModel<<<NBLOCKS,NTHREADS >>>(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, 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){ 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<<<NBLOCKS,NTHREADS >>>(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, 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){ 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<<<NBLOCKS,NTHREADS >>>(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){
}

126
models/ColorModel.cpp
View File

@ -531,6 +531,121 @@ void ScaLBL_ColorModel::Initialize(){
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); 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<Database> 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<double>( 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(){ void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz; int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,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" )){ if (color_db->keyExists( "krA_morph_factor" )){
KRA_MORPH_FACTOR = color_db->getScalar<double>( "krA_morph_factor" ); KRA_MORPH_FACTOR = color_db->getScalar<double>( "krA_morph_factor" );
} }
/* defaults for simulation protocols */ /* defaults for simulation protocols */
auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" ); auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){ if (protocol == "image sequence"){
@ -625,7 +739,7 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "seed_water" )){ if (analysis_db->keyExists( "seed_water" )){
seed_water = analysis_db->getScalar<double>( "seed_water" ); seed_water = analysis_db->getScalar<double>( "seed_water" );
if (rank == 0) printf("Seed water in oil %f (seed_water) \n",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" )){ if (analysis_db->keyExists( "morph_delta" )){
morph_delta = analysis_db->getScalar<double>( "morph_delta" ); morph_delta = analysis_db->getScalar<double>( "morph_delta" );
@ -656,7 +770,6 @@ void ScaLBL_ColorModel::Run(){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" ); MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
} }
if (rank==0){ if (rank==0){
printf("********************************************************\n"); printf("********************************************************\n");
if (protocol == "image sequence"){ if (protocol == "image sequence"){
@ -1320,7 +1433,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
double vF = 0.f; double vF = 0.f;
double vS = 0.f; double vS = 0.f;
double delta_volume; double delta_volume;
double WallFactor = 0.0; double WallFactor = 1.0;
bool USE_CONNECTED_NWP = false; bool USE_CONNECTED_NWP = false;
DoubleArray phase(Nx,Ny,Nz); 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 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); sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
FILE *INPUT = fopen(LocalRankFilename,"wb"); FILE *INPUT = fopen(LocalRankFilename,"wb");

6
models/ColorModel.h
View File

@ -16,6 +16,10 @@ Implementation of color lattice boltzmann model
#include "ProfilerApp.h" #include "ProfilerApp.h"
#include "threadpool/thread_pool.h" #include "threadpool/thread_pool.h"
#ifndef ScaLBL_ColorModel_INC
#define ScaLBL_ColorModel_INC
class ScaLBL_ColorModel{ class ScaLBL_ColorModel{
public: public:
ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM); ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM);
@ -29,6 +33,7 @@ public:
void Create(); void Create();
void Initialize(); void Initialize();
void Run(); void Run();
double Run(int returntime);
void WriteDebug(); void WriteDebug();
void getPhaseField(DoubleArray &f); void getPhaseField(DoubleArray &f);
@ -99,4 +104,5 @@ private:
int timestep; int timestep;
int timestep_previous; int timestep_previous;
}; };
#endif

333
models/FreeLeeModel.cpp
View File

@ -10,9 +10,9 @@ color lattice boltzmann model
#include <time.h> #include <time.h>
ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM): 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), 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) 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(){ 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; k<Nzh-1; k++){
for (int j=1; j<Nyh-1; j++){
for (int i=1; i<Nxh-1; i++){
PhaseValues(i-1,j-1,k-1) = PhaseWideHalo(i,j,k);
}
}
}
}
void ScaLBL_FreeLeeModel::getPotential(DoubleArray &PressureValues, DoubleArray &MuValues){
ScaLBL_Comm->RegularLayout(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){ void ScaLBL_FreeLeeModel::ReadParams(string filename){
// read the input database // read the input database
db = std::make_shared<Database>( filename ); db = std::make_shared<Database>( filename );
@ -529,6 +568,34 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
//fwrite(phase,8,Nh,OUTFILE); //fwrite(phase,8,Nh,OUTFILE);
//fclose(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 [] phase;
delete [] ColorGrad_host; 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; 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 ***************************************/ //************ MAIN ITERATION LOOP ***************************************/
comm.barrier(); comm.barrier();
auto t1 = std::chrono::system_clock::now(); auto t1 = std::chrono::system_clock::now();
PROFILE_START("Loop"); 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)); } //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update"); PROFILE_START("Update");
// *************ODD TIMESTEP************* // *************ODD TIMESTEP*************
@ -732,24 +795,27 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){
// Compute the Phase indicator field // Compute the Phase indicator field
// Read for hq happens in this routine (requires communication) // Read for hq happens in this routine (requires communication)
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL 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->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); 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 // 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){ if (BoundaryCondition > 0 && BoundaryCondition < 5){
//TODO to be revised //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, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM, ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
ScaLBL_Comm_WideHalo->Recv(Phi); 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->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
// Set BCs // 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_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_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(); ScaLBL_Comm->Barrier();
// *************EVEN TIMESTEP************* // *************EVEN TIMESTEP*************
timestep++; timestep++;
// Compute the Phase indicator field // Compute the Phase indicator field
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMA
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); 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->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); 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 // Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
// Halo exchange for phase field // 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){ if (BoundaryCondition > 0 && BoundaryCondition < 5){
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
ScaLBL_Comm_WideHalo->Send(Phi); ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
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_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
ScaLBL_Comm_WideHalo->Recv(Phi); 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->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
// Set boundary conditions // 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_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, ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
//************************************************************************ //************************************************************************
PROFILE_STOP("Update"); PROFILE_STOP("Update");
@ -816,18 +886,11 @@ void ScaLBL_FreeLeeModel::Run_TwoFluid(){
if (rank==0) printf("-------------------------------------------------------------------\n"); if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep // Compute the walltime per timestep
auto t2 = std::chrono::system_clock::now(); auto t2 = std::chrono::system_clock::now();
double cputime = std::chrono::duration<double>( t2 - t1 ).count() / timestep; double cputime = std::chrono::duration<double>( t2 - t1 ).count() / (EXIT_TIME-START_TIME);
// Performance obtained from each node // Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000; double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n"); return MLUPS;
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");
// ************************************************************************
} }
void ScaLBL_FreeLeeModel::Run_SingleFluid(){ void ScaLBL_FreeLeeModel::Run_SingleFluid(){
@ -878,6 +941,7 @@ void ScaLBL_FreeLeeModel::Run_SingleFluid(){
0, ScaLBL_Comm->LastExterior(), Np); 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
// *************EVEN TIMESTEP************* // *************EVEN TIMESTEP*************
timestep++; timestep++;
//------------------------------------------------------------------------------------------------------------------- //-------------------------------------------------------------------------------------------------------------------
@ -932,6 +996,32 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
DoubleArray PhaseData(Nxh,Nyh,Nzh); DoubleArray PhaseData(Nxh,Nyh,Nzh);
//ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
ScaLBL_CopyToHost(PhaseData.data(), Phi, sizeof(double)*Nh); 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; FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.%05i.raw",rank); sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
@ -940,6 +1030,17 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
fclose(OUTFILE); fclose(OUTFILE);
DoubleArray PhaseField(Nx,Ny,Nz); 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); ScaLBL_Comm->RegularLayout(Map,Den,PhaseField);
FILE *AFILE; FILE *AFILE;
sprintf(LocalRankFilename,"Density.%05i.raw",rank); sprintf(LocalRankFilename,"Density.%05i.raw",rank);
@ -975,7 +1076,7 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
fwrite(PhaseField.data(),8,N,VELZ_FILE); fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE); fclose(VELZ_FILE);
/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField); ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
FILE *CGX_FILE; FILE *CGX_FILE;
sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank); sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
CGX_FILE = fopen(LocalRankFilename,"wb"); CGX_FILE = fopen(LocalRankFilename,"wb");
@ -995,7 +1096,7 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
CGZ_FILE = fopen(LocalRankFilename,"wb"); CGZ_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGZ_FILE); fwrite(PhaseField.data(),8,N,CGZ_FILE);
fclose(CGZ_FILE); fclose(CGZ_FILE);
*/
} }
void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){ void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){
@ -1031,3 +1132,151 @@ void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){
fwrite(PhaseField.data(),8,N,VELZ_FILE); fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE); fclose(VELZ_FILE);
} }
void ScaLBL_FreeLeeModel::Create_DummyPhase_MGTest(){
// Initialize communication structures in averaging domain
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[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<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
//ScaLBL_Comm_Regular = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
ScaLBL_Comm_WideHalo = std::shared_ptr<ScaLBLWideHalo_Communicator>(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; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
TmpMap[idx] = ScaLBL_Comm_WideHalo->Map(i,j,k);
}
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); 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(); idx<ScaLBL_Comm->LastInterior(); 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;k<Nzh;k++){
for (int j=0;j<Nyh;j++){
for (int i=0;i<Nxh;i++){
//idx for double-halo array 'phase'
int nh = k*Nxh*Nyh+j*Nxh+i;
//idx for single-halo array Mask->id[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(); i<ScaLBL_Comm->LastInterior();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; i<ScaLBL_Comm->LastExterior();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;
}

16
models/FreeLeeModel.h
View File

@ -16,6 +16,9 @@ Implementation of Lee et al JCP 2016 lattice boltzmann model
#include "common/ScaLBL.h" #include "common/ScaLBL.h"
#include "common/WideHalo.h" #include "common/WideHalo.h"
#ifndef ScaLBL_FreeLeeModel_INC
#define ScaLBL_FreeLeeModel_INC
class ScaLBL_FreeLeeModel{ class ScaLBL_FreeLeeModel{
public: public:
ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM); ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM);
@ -28,12 +31,17 @@ public:
void ReadInput(); void ReadInput();
void Create_TwoFluid(); void Create_TwoFluid();
void Initialize_TwoFluid(); void Initialize_TwoFluid();
void Run_TwoFluid(); double Run_TwoFluid(int returntime);
void WriteDebug_TwoFluid(); void WriteDebug_TwoFluid();
void Create_SingleFluid(); void Create_SingleFluid();
void Initialize_SingleFluid(); void Initialize_SingleFluid();
void Run_SingleFluid(); void Run_SingleFluid();
void WriteDebug_SingleFluid(); void WriteDebug_SingleFluid();
// test utilities
void Create_DummyPhase_MGTest();
void MGTest();
bool Restart,pBC; bool Restart,pBC;
int timestep,timestepMax; int timestep,timestepMax;
@ -73,6 +81,10 @@ public:
double *Velocity; double *Velocity;
double *Pressure; double *Pressure;
void getPhase(DoubleArray &PhaseValues);
void getPotential(DoubleArray &PressureValues, DoubleArray &MuValues);
void getVelocity(DoubleArray &Vx, DoubleArray &Vy, DoubleArray &Vz);
DoubleArray SignDist; DoubleArray SignDist;
private: private:
@ -90,4 +102,4 @@ private:
void AssignComponentLabels_ChemPotential_ColorGrad(); void AssignComponentLabels_ChemPotential_ColorGrad();
}; };
#endif

3
sample_scripts/configure_ubuntu
View File

@ -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_FLAGS="-O3 -fPIC " \
-D CMAKE_CXX_STANDARD=14 \ -D CMAKE_CXX_STANDARD=14 \
-D MPIEXEC=mpirun \ -D MPIEXEC=mpirun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D CMAKE_BUILD_TYPE:STRING=Release \ -D CMAKE_BUILD_TYPE:STRING=Release \
-D CUDA_FLAGS="-arch sm_35" \ -D CUDA_FLAGS="-arch sm_35" \
-D CUDA_HOST_COMPILER="/usr/bin/gcc" \ -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 USE_SILO=1 \
-D SILO_LIB="/opt/arden/silo/4.10.2/lib/libsiloh5.a" \ -D SILO_LIB="/opt/arden/silo/4.10.2/lib/libsiloh5.a" \
-D SILO_DIRECTORY="/opt/arden/silo/4.10.2" \ -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 NETCDF_DIRECTORY="/opt/arden/netcdf/4.6.1" \
-D USE_CUDA=0 \ -D USE_CUDA=0 \
-D USE_TIMER=0 \ -D USE_TIMER=0 \

2
tests/CMakeLists.txt
View File

@ -43,6 +43,7 @@ ADD_LBPM_EXECUTABLE( TestPoissonSolver )
ADD_LBPM_EXECUTABLE( TestIonModel ) ADD_LBPM_EXECUTABLE( TestIonModel )
ADD_LBPM_EXECUTABLE( TestNernstPlanck ) ADD_LBPM_EXECUTABLE( TestNernstPlanck )
ADD_LBPM_EXECUTABLE( TestPNP_Stokes ) ADD_LBPM_EXECUTABLE( TestPNP_Stokes )
ADD_LBPM_EXECUTABLE( TestMixedGrad )
@ -61,6 +62,7 @@ ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestWideHalo ) ADD_LBPM_TEST( TestWideHalo )
ADD_LBPM_TEST( TestColorGradDFH ) ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db) 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( TestColorMassBounceback ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db) ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db) ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db)

199
tests/TestMixedGrad.cpp Normal file
View File

@ -0,0 +1,199 @@
#include <exception>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#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;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
LeeModel.Mask->id[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;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
LeeModel.Mask->id[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(); n<LeeModel.ScaLBL_Comm->LastInterior(); 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<Database>( 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;
}

34
tests/lbpm_color_simulator.cpp
View File

@ -27,19 +27,24 @@ int main( int argc, char **argv )
// Initialize // Initialize
Utilities::startup( argc, argv ); Utilities::startup( argc, argv );
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
{ // Limit scope so variables that contain communicators will free before MPI_Finialize { // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::MPI comm( MPI_COMM_WORLD ); Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank(); int rank = comm.getRank();
int nprocs = comm.getSize(); int nprocs = comm.getSize();
std::string SimulationMode = "production";
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
if (argc > 2) {
SimulationMode = "development";
}
if ( rank == 0 ) { if ( rank == 0 ) {
printf( "********************************************************\n" ); printf( "********************************************************\n" );
printf( "Running Color LBM \n" ); printf( "Running Color LBM \n" );
printf( "********************************************************\n" ); printf( "********************************************************\n" );
if (SimulationMode == "development")
printf("**** DEVELOPMENT MODE ENABLED *************\n");
} }
// Initialize compute device // Initialize compute device
int device = ScaLBL_SetDevice( rank ); 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 ColorModel.Create(); // creating the model will create data structure to match the pore
// structure and allocate variables // structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables
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<int>( "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.Run();
// ColorModel.WriteDebug();
ColorModel.WriteDebug();
PROFILE_STOP( "Main" ); PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" ); auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );

28
tests/lbpm_freelee_simulator.cpp
View File

@ -8,6 +8,7 @@
#include "common/Utilities.h" #include "common/Utilities.h"
#include "models/FreeLeeModel.h" #include "models/FreeLeeModel.h"
#include "analysis/FreeEnergy.h"
//******************************************************************* //*******************************************************************
// Implementation of Free-Energy Two-Phase LBM (Lee model) // Implementation of Free-Energy Two-Phase LBM (Lee model)
@ -52,9 +53,32 @@ int main( int argc, char **argv )
LeeModel.SetDomain(); LeeModel.SetDomain();
LeeModel.ReadInput(); LeeModel.ReadInput();
LeeModel.Create_TwoFluid(); LeeModel.Create_TwoFluid();
FreeEnergyAnalyzer Analysis(LeeModel.Dm);
LeeModel.Initialize_TwoFluid(); 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<int>( "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"); PROFILE_STOP("Main");
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_simulator" ); auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_simulator" );

101
tests/testGlobalMassFreeLee.cpp Normal file
View File

@ -0,0 +1,101 @@
#include <exception>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#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<Database>( 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<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
double change = DensityFinal(i,j,k)-DensityInit(i,j,k);
DensityChange(i,j,k) = change;
totalChange += change;
}
}
}
printf("Density change, %f\n", totalChange);
FILE *OUTFILE;
char LocalRankFilename[40];
sprintf(LocalRankFilename,"DensityChange.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(DensityChange.data(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
//LeeModel.WriteDebug_TwoFluid();
PROFILE_STOP("Main");
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
return 0;
}