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

Merge branch 'tmpfix' into FOM

Browse Source
This commit is contained in:
James McClure
2021-01-15 17:13:08 -05:00
parent 31234ed003 c911a4f3d2
commit 2af6aff545
29 changed files with 4276 additions and 243 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
common/Domain.cpp
View File

@@ -1,5 +1,5 @@
// Created by James McClure
// Copyright 2008-2013
// Copyright 2008-2020
#include <stdio.h>
#include <stdlib.h>
#include <iostream>

23
common/ScaLBL.cpp
View File

@@ -351,8 +351,7 @@ void ScaLBL_Communicator::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *li
delete [] ReturnDist;
}
int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np){
int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width){
/*
* Generate a memory optimized layout
* id[n] == 0 implies that site n should be ignored (treat as a mask)
@@ -395,28 +394,26 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
n = k*Nx*Ny+j*Nx+i;
if (id[n] > 0){
// Counts for the six faces
if (i==1) Map(n)=idx++;
else if (j==1) Map(n)=idx++;
else if (k==1) Map(n)=idx++;
else if (i==Nx-2) Map(n)=idx++;
else if (j==Ny-2) Map(n)=idx++;
else if (k==Nz-2) Map(n)=idx++;
if (i>0 && i<=width) Map(n)=idx++;
else if (j>0 && j<=width) Map(n)=idx++;
else if (k>0 && k<=width) Map(n)=idx++;
else if (i>Nx-width-2 && i<Nx-1) Map(n)=idx++;
else if (j>Ny-width-2 && j<Ny-1) Map(n)=idx++;
else if (k>Nz-width-2 && k<Nz-1) Map(n)=idx++;
}
}
}
}
next=idx;
//printf("Interior... \n");
// ********* Interior **********
// align the next read
first_interior=(next/16 + 1)*16;
idx = first_interior;
// Step 2/2: Next loop over the domain interior in block-cyclic fashion
for (k=2; k<Nz-2; k++){
for (j=2; j<Ny-2; j++){
for (i=2; i<Nx-2; i++){
for (k=width+1; k<Nz-width-1; k++){
for (j=width+1; j<Ny-width-1; j++){
for (i=width+1; i<Nx-width-1; i++){
// Local index (regular layout)
n = k*Nx*Ny + j*Nx + i;
if (id[n] > 0 ){

145
common/ScaLBL.h
View File

@@ -119,104 +119,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int f
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,
double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
// GREYSCALE FREE-ENERGY MODEL (Two-component)
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleFE(double *dist, double *Aq, double *Bq, double *Den,
// double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
// double *Poros,double *Perm, double *Velocity,double *Pressure);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleFE(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,
// double *DenGradA, double *DenGradB, double *SolidForce, int start, int finish, int Np,
// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz,
// double *Poros,double *Perm, double *Velocity,double *Pressure);
//
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleFEChem(double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np,
// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB,
// double Gx, double Gy, double Gz,
// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleFEChem(int *neighborList, double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np,
// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB,
// double Gx, double Gy, double Gz,
// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap);
//
//extern "C" void ScaLBL_D3Q7_GreyscaleFE_Init(double *Den, double *Cq, double *PhiLap, double gamma, double kappaA, double kappaB, double lambdaA, double lambdaB, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleFE_IMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np);
//
//extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleFEDensity(int *NeighborList, double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleFEDensity(double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleFEPhi(int *NeighborList, double *Cq, double *Phi, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleFEPhi(double *Cq, double *Phi, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Laplacian(int *neighborList, double *Den, double *DenLap, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleFE_Pressure(double *dist, double *Den, double *Porosity,double *Velocity,
// double *Pressure, double rhoA,double rhoB, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleFE_PressureTensor(int *neighborList, double *Phi,double *Pressure, double *PressTensor, double *PhiLap,
// double kappaA,double kappaB,double lambdaA,double lambdaB, int start, int finish, int Np);
// GREYSCALE SHAN-CHEN MODEL (Two-component)
//extern "C" void ScaLBL_D3Q19_GreyscaleSC_Init(int *Map, double *distA, double *distB, double *DenA, double *DenB, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, int *Mpa, double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB,
// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure,
// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz,
// int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(int *Map,double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB,
// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure,
// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz,
// int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB,
// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure,
// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz,
// int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB,
// double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure,
// double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz,
// int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_GreyscaleSC_Gradient(int *neighborList, int *Map, double *Den, double *DenGrad, int strideY, int strideZ,int start, int finish, int Np);
//
//extern "C" void ScaLBL_GreyscaleSC_BC_z(int *list, int *Map, double *DenA, double *DenB, double vA, double vB, int count);
//
//extern "C" void ScaLBL_GreyscaleSC_BC_Z(int *list, int *Map, double *DenA, double *DenB, double vA, double vB, int count);
//
//extern "C" void ScaLBL_GreyscaleSC_AAeven_Pressure_BC_z(int *list, double *distA, double *distB, double dinA, double dinB, int count, int N);
//
//extern "C" void ScaLBL_GreyscaleSC_AAeven_Pressure_BC_Z(int *list, double *distA, double *distB, double doutA, double doutB, int count, int N);
//
//extern "C" void ScaLBL_GreyscaleSC_AAodd_Pressure_BC_z(int *neighborList, int *list, double *distA, double *distB, double dinA, double dinB, int count, int N);
//
//extern "C" void ScaLBL_GreyscaleSC_AAodd_Pressure_BC_Z(int *neighborList, int *list, double *distA, double *distB, double doutA, double doutB, int count, int N);
// GREYSCALE COLOR MODEL (Two-component)
//extern "C" void ScaLBL_D3Q19_GreyscaleColor_Init(double *dist, double *Porosity, int Np);
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
// double *ColorGrad,double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
//
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
// double *ColorGrad,double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
// double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure,
@@ -227,6 +129,49 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
// ION TRANSPORT MODEL
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np);
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);
// LBM Poisson solver
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
//maybe deprecated
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np);
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);
@@ -254,6 +199,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq,
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_PhaseField_Init(int *Map, double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np);
// Density functional hydrodynamics LBM
@@ -370,7 +317,7 @@ public:
int FirstInterior();
int LastInterior();
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np);
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width);
void Barrier(){
ScaLBL_DeviceBarrier();
MPI_COMM_SCALBL.barrier();

338
common/WideHalo.cpp Normal file
View File

@@ -0,0 +1,338 @@
/*
This class implements support for halo widths larger than 1
*/
#include "common/WideHalo.h"
ScaLBLWideHalo_Communicator::ScaLBLWideHalo_Communicator(std::shared_ptr <Domain> Dm, int width)
{
//......................................................................................
Lock=false; // unlock the communicator
//......................................................................................
// Create a separate copy of the communicator for the device
MPI_COMM_SCALBL = Dm->Comm.dup();
//......................................................................................
// Copy the domain size and communication information directly from Dm
Nx = Dm->Nx;
Ny = Dm->Ny;
Nz = Dm->Nz;
N = Nx*Ny*Nz;
Nxh = Nx + 2*(width - 1);
Nyh = Ny + 2*(width - 1);
Nzh = Nz + 2*(width - 1);
Nh = Nxh*Nyh*Nzh;
Map.resize(Nx,Ny,Nz);
rank=Dm->rank();
iproc = Dm->iproc();
jproc = Dm->jproc();
kproc = Dm->kproc();
nprocx = Dm->nprocx();
nprocy = Dm->nprocy();
nprocz = Dm->nprocz();
rank_info = RankInfoStruct(rank,nprocx,nprocy,nprocz);
rank = rank_info.rank[1][1][1];
rank_X = rank_info.rank[2][1][1];
rank_x = rank_info.rank[0][1][1];
rank_Y = rank_info.rank[1][2][1];
rank_y = rank_info.rank[1][0][1];
rank_Z = rank_info.rank[1][1][2];
rank_z = rank_info.rank[1][1][0];
rank_XY = rank_info.rank[2][2][1];
rank_xy = rank_info.rank[0][0][1];
rank_Xy = rank_info.rank[2][0][1];
rank_xY = rank_info.rank[0][2][1];
rank_XZ = rank_info.rank[2][1][2];
rank_xz = rank_info.rank[0][1][0];
rank_Xz = rank_info.rank[2][1][0];
rank_xZ = rank_info.rank[0][1][2];
rank_YZ = rank_info.rank[1][2][2];
rank_yz = rank_info.rank[1][0][0];
rank_Yz = rank_info.rank[1][2][0];
rank_yZ = rank_info.rank[1][0][2];
rank_XYz = rank_info.rank[2][2][0];
rank_xyz = rank_info.rank[0][0][0];
rank_Xyz = rank_info.rank[2][0][0];
rank_xYz = rank_info.rank[0][2][0];
rank_XYZ = rank_info.rank[2][2][2];
rank_xyZ = rank_info.rank[0][0][2];
rank_XyZ = rank_info.rank[2][0][2];
rank_xYZ = rank_info.rank[0][2][2];
MPI_COMM_SCALBL.barrier();
/* Fill in communications patterns for the lists */
/* Send lists */
sendCount_x = getHaloBlock(width,2*width,width,Nyh-width,width,Nzh-width,dvcSendList_x);
sendCount_X =getHaloBlock(Nxh-2*width,Nxh-width,width,Nyh-width,width,Nzh-width,dvcSendList_X);
sendCount_y =getHaloBlock(width,Nxh-width,width,2*width,width,Nzh-width,dvcSendList_y);
sendCount_Y =getHaloBlock(width,Nxh-width,Nyh-2*width,Nyh-width,width,Nzh-width,dvcSendList_Y);
sendCount_z =getHaloBlock(width,Nxh-width,width,Nyh-width,width,2*width,dvcSendList_z);
sendCount_X =getHaloBlock(width,Nxh-width,width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_Z);
// xy
sendCount_xy =getHaloBlock(width,2*width,width,2*width,width,Nzh-width,dvcSendList_xy);
sendCount_xY =getHaloBlock(width,2*width,Nyh-2*width,Nyh-width,width,Nzh-width,dvcSendList_xY);
sendCount_Xy =getHaloBlock(Nxh-2*width,Nxh-width,width,2*width,width,Nzh-width,dvcSendList_Xy);
sendCount_XY =getHaloBlock(Nxh-2*width,Nxh-width,Nyh-2*width,Nyh-width,width,Nzh-width,dvcSendList_XY);
// xz
sendCount_xz =getHaloBlock(width,2*width,width,Nyh-width,width,2*width,dvcSendList_xz);
sendCount_xZ =getHaloBlock(width,2*width,width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_xZ);
sendCount_Xz =getHaloBlock(Nxh-2*width,Nxh-width,width,Nyh-width,width,2*width,dvcSendList_Xz);
sendCount_XZ =getHaloBlock(Nxh-2*width,Nxh-width,width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_XZ);
// yz
sendCount_yz =getHaloBlock(width,Nxh-width,width,2*width,width,2*width,dvcSendList_yz);
sendCount_yZ =getHaloBlock(width,Nxh-width,width,2*width,Nzh-2*width,Nzh-width,dvcSendList_yZ);
sendCount_Yz =getHaloBlock(width,Nxh-width,Nyh-2*width,Nyh-width,width,2*width,dvcSendList_Yz);
sendCount_YZ =getHaloBlock(width,Nxh-width,Nyh-2*width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_YZ);
// xyz
sendCount_xyz =getHaloBlock(width,2*width,width,2*width,width,2*width,dvcSendList_xyz);
sendCount_xyZ =getHaloBlock(width,2*width,width,2*width,Nzh-2*width,Nzh-width,dvcSendList_xyZ);
sendCount_xYz =getHaloBlock(width,2*width,Nyh-2*width,Nyh-width,width,2*width,dvcSendList_xYz);
sendCount_xYz =getHaloBlock(width,2*width,Nyh-2*width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_xYZ);
sendCount_Xyz =getHaloBlock(Nxh-2*width,Nxh-width,width,2*width,width,2*width,dvcSendList_Xyz);
sendCount_XyZ =getHaloBlock(Nxh-2*width,Nxh-width,width,2*width,Nzh-2*width,Nzh-width,dvcSendList_XyZ);
sendCount_XYz =getHaloBlock(Nxh-2*width,Nxh-width,Nyh-2*width,Nyh-width,width,2*width,dvcSendList_XYz);
sendCount_XYZ =getHaloBlock(Nxh-2*width,Nxh-width,Nyh-2*width,Nyh-width,Nzh-2*width,Nzh-width,dvcSendList_XYZ);
/* Recv lists */
recvCount_x =getHaloBlock(0,width,width,Nyh-width,width,Nzh-width,dvcRecvList_x);
recvCount_X =getHaloBlock(Nxh-width,Nxh,width,Nyh-width,width,Nzh-width,dvcRecvList_X);
recvCount_y =getHaloBlock(width,Nxh-width,0,width,width,Nzh-width,dvcRecvList_y);
recvCount_Y =getHaloBlock(width,Nxh-width,Nyh-width,Nyh,width,Nzh-width,dvcRecvList_Y);
recvCount_z =getHaloBlock(width,Nxh-width,width,Nyh-width,0,width,dvcRecvList_z);
recvCount_Z =getHaloBlock(width,Nxh-width,width,Nyh-width,Nzh-width,Nzh,dvcRecvList_Z);
//xy
recvCount_xy =getHaloBlock(0,width,0,width,width,Nzh-width,dvcRecvList_xy);
recvCount_xY =getHaloBlock(0,width,Nyh-width,Nyh,width,Nzh-width,dvcRecvList_xY);
recvCount_Xy =getHaloBlock(Nxh-width,Nxh,0,width,width,Nzh-width,dvcRecvList_Xy);
recvCount_XY =getHaloBlock(Nxh-width,Nxh,Nyh-width,Nyh,width,Nzh-width,dvcRecvList_XY);
//xz
recvCount_xz =getHaloBlock(0,width,width,Nyh-width,0,width,dvcRecvList_xz);
recvCount_xZ =getHaloBlock(0,width,width,Nyh-width,Nzh-width,Nzh,dvcRecvList_xZ);
recvCount_Xz =getHaloBlock(Nxh-width,Nxh,width,Nyh-width,0,width,dvcRecvList_Xz);
recvCount_XZ =getHaloBlock(Nxh-width,Nxh,width,Nyh-width,Nzh-width,Nzh,dvcRecvList_XZ);
//yz
recvCount_yz =getHaloBlock(width,Nxh-width,0,width,0,width,dvcRecvList_yz);
recvCount_yZ =getHaloBlock(width,Nxh-width,0,width,Nzh-width,Nzh,dvcRecvList_yZ);
recvCount_Yz =getHaloBlock(width,Nxh-width,Nyh-width,Nyh,0,width,dvcRecvList_Yz);
recvCount_YZ =getHaloBlock(width,Nxh-width,Nyh-width,Nyh,Nzh-width,Nzh,dvcRecvList_YZ);
//xyz
recvCount_xyz =getHaloBlock(0,width,0,width,0,width,dvcRecvList_xyz);
recvCount_xyZ =getHaloBlock(0,width,0,width,Nzh-width,Nzh,dvcRecvList_xyZ);
recvCount_xYz =getHaloBlock(0,width,Nyh-width,Nyh,0,width,dvcRecvList_xYz);
recvCount_xYZ =getHaloBlock(0,width,Nyh-width,Nyh,Nzh-width,Nzh,dvcRecvList_xYZ);
recvCount_Xyz =getHaloBlock(Nxh-width,Nxh,0,width,0,width,dvcRecvList_Xyz);
recvCount_XyZ =getHaloBlock(Nxh-width,Nxh,0,width,Nzh-width,Nzh,dvcRecvList_XyZ);
recvCount_XYz =getHaloBlock(Nxh-width,Nxh,Nyh-width,Nyh,0,width,dvcRecvList_XYz);
recvCount_XYZ =getHaloBlock(Nxh-width,Nxh,Nyh-width,Nyh,Nzh-width,Nzh,dvcRecvList_XYZ);
//......................................................................................
ScaLBL_AllocateZeroCopy((void **) &sendbuf_x, sendCount_x*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_X, sendCount_X*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_y, sendCount_y*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Y, sendCount_Y*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_z, sendCount_z*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Z, sendCount_Z*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xy, sendCount_xy*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xY, sendCount_xY*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xy, sendCount_Xy*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_XY, sendCount_XY*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xz, sendCount_xz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xZ, sendCount_xZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xz, sendCount_Xz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_XZ, sendCount_XZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_yz, sendCount_yz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_yZ, sendCount_yZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Yz, sendCount_Yz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_YZ, sendCount_YZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xyz, sendCount_xyz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xYz, sendCount_xYz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_Xyz, sendCount_Xyz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_XYz, sendCount_XYz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xyZ, sendCount_xyZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_xYZ, sendCount_xYZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_XyZ, sendCount_XyZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &sendbuf_XYZ, sendCount_XYZ*sizeof(double)); // Allocate device memory
//......................................................................................
ScaLBL_AllocateZeroCopy((void **) &recvbuf_x, recvCount_x*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_X, recvCount_X*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_y, recvCount_y*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Y, recvCount_Y*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_z, recvCount_z*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Z, recvCount_Z*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xy, recvCount_xy*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xY, recvCount_xY*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xy, recvCount_Xy*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_XY, recvCount_XY*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xz, recvCount_xz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xZ, recvCount_xZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xz, recvCount_Xz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_XZ, recvCount_XZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_yz, recvCount_yz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_yZ, recvCount_yZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Yz, recvCount_Yz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_YZ, recvCount_YZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xyz, recvCount_xyz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xYz, recvCount_xYz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_Xyz, recvCount_Xyz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_XYz, recvCount_XYz*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xyZ, recvCount_xyZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_xYZ, recvCount_xYZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_XyZ, recvCount_XyZ*sizeof(double)); // Allocate device memory
ScaLBL_AllocateZeroCopy((void **) &recvbuf_XYZ, recvCount_XYZ*sizeof(double)); // Allocate device memory
/* Set up a map to the halo width=1 data structure */
for (k=width; k<Nzh-width; k++){
for (j=width; j<Nyh-width; j++){
for (i=width; i<Nxh-width; i++){
int idx = k*Nxh*Nyh + j*Nxh + i;
Map(i-width+1,j-width+1,k-width+1) = idx;
}
}
}
}
void ScaLBLWideHalo_Communicator::Send(double *data){
//...................................................................................
if (Lock==true){
ERROR("ScaLBL Error (SendHalo): ScaLBLWideHalo_Communicator is locked -- did you forget to match Send/Recv calls?");
}
else{
Lock=true;
}
ScaLBL_DeviceBarrier();
//...................................................................................
sendtag = recvtag = 1;
//...................................................................................
ScaLBL_Scalar_Pack(dvcSendList_x, sendCount_x,sendbuf_x, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_y, sendCount_y,sendbuf_y, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_z, sendCount_z,sendbuf_z, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_X, sendCount_X,sendbuf_X, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Y, sendCount_Y,sendbuf_Y, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Z, sendCount_Z,sendbuf_Z, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xy, sendCount_xy,sendbuf_xy, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xY, sendCount_xY,sendbuf_xY, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Xy, sendCount_Xy,sendbuf_Xy, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_XY, sendCount_XY,sendbuf_XY, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xz, sendCount_xz,sendbuf_xz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xZ, sendCount_xZ,sendbuf_xZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Xz, sendCount_Xz,sendbuf_Xz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_XZ, sendCount_XZ,sendbuf_XZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_yz, sendCount_yz,sendbuf_yz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_yZ, sendCount_yZ,sendbuf_yZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Yz, sendCount_Yz,sendbuf_Yz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_YZ, sendCount_YZ,sendbuf_YZ, data, Nh);
/* corners */
ScaLBL_Scalar_Pack(dvcSendList_xyz, sendCount_xyz,sendbuf_xyz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xyZ, sendCount_xyZ,sendbuf_xyZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xYz, sendCount_xYz,sendbuf_xYz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_xYZ, sendCount_xYZ,sendbuf_xYZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_Xyz, sendCount_Xyz,sendbuf_Xyz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_XyZ, sendCount_XyZ,sendbuf_XyZ, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_XYz, sendCount_XYz,sendbuf_XYz, data, Nh);
ScaLBL_Scalar_Pack(dvcSendList_XYZ, sendCount_XYZ,sendbuf_XYZ, data, Nh);
//...................................................................................
// Send / Recv all the phase indcator field values
//...................................................................................
req1[0] = MPI_COMM_SCALBL.Isend(&sendCount_x,1,rank_x,sendtag+0);
req2[0] = MPI_COMM_SCALBL.Irecv(&recvCount_X,1,rank_X,recvtag+0);
req1[1] = MPI_COMM_SCALBL.Isend(&sendCount_X,1,rank_X,sendtag+1);
req2[1] = MPI_COMM_SCALBL.Irecv(&recvCount_x,1,rank_x,recvtag+1);
req1[2] = MPI_COMM_SCALBL.Isend(&sendCount_y,1,rank_y,sendtag+2);
req2[2] = MPI_COMM_SCALBL.Irecv(&recvCount_Y,1,rank_Y,recvtag+2);
req1[3] = MPI_COMM_SCALBL.Isend(&sendCount_Y,1,rank_Y,sendtag+3);
req2[3] = MPI_COMM_SCALBL.Irecv(&recvCount_y,1,rank_y,recvtag+3);
req1[4] = MPI_COMM_SCALBL.Isend(&sendCount_z,1,rank_z,sendtag+4);
req2[4] = MPI_COMM_SCALBL.Irecv(&recvCount_Z,1,rank_Z,recvtag+4);
req1[5] = MPI_COMM_SCALBL.Isend(&sendCount_Z,1,rank_Z,sendtag+5);
req2[5] = MPI_COMM_SCALBL.Irecv(&recvCount_z,1,rank_z,recvtag+5);
req1[6] = MPI_COMM_SCALBL.Isend(&sendCount_xy,1,rank_xy,sendtag+6);
req2[6] = MPI_COMM_SCALBL.Irecv(&recvCount_XY,1,rank_XY,recvtag+6);
req1[7] = MPI_COMM_SCALBL.Isend(&sendCount_XY,1,rank_XY,sendtag+7);
req2[7] = MPI_COMM_SCALBL.Irecv(&recvCount_xy,1,rank_xy,recvtag+7);
req1[8] = MPI_COMM_SCALBL.Isend(&sendCount_Xy,1,rank_Xy,sendtag+8);
req2[8] = MPI_COMM_SCALBL.Irecv(&recvCount_xY,1,rank_xY,recvtag+8);
req1[9] = MPI_COMM_SCALBL.Isend(&sendCount_xY,1,rank_xY,sendtag+9);
req2[9] = MPI_COMM_SCALBL.Irecv(&recvCount_Xy,1,rank_Xy,recvtag+9);
req1[10] = MPI_COMM_SCALBL.Isend(&sendCount_xz,1,rank_xz,sendtag+10);
req2[10] = MPI_COMM_SCALBL.Irecv(&recvCount_XZ,1,rank_XZ,recvtag+10);
req1[11] = MPI_COMM_SCALBL.Isend(&sendCount_XZ,1,rank_XZ,sendtag+11);
req2[11] = MPI_COMM_SCALBL.Irecv(&recvCount_xz,1,rank_xz,recvtag+11);
req1[12] = MPI_COMM_SCALBL.Isend(&sendCount_Xz,1,rank_Xz,sendtag+12);
req2[12] = MPI_COMM_SCALBL.Irecv(&recvCount_xZ,1,rank_xZ,recvtag+12);
req1[13] = MPI_COMM_SCALBL.Isend(&sendCount_xZ,1,rank_xZ,sendtag+13);
req2[13] = MPI_COMM_SCALBL.Irecv(&recvCount_Xz,1,rank_Xz,recvtag+13);
req1[14] = MPI_COMM_SCALBL.Isend(&sendCount_yz,1,rank_yz,sendtag+14);
req2[14] = MPI_COMM_SCALBL.Irecv(&recvCount_YZ,1,rank_YZ,recvtag+14);
req1[15] = MPI_COMM_SCALBL.Isend(&sendCount_YZ,1,rank_YZ,sendtag+15);
req2[15] = MPI_COMM_SCALBL.Irecv(&recvCount_yz,1,rank_yz,recvtag+15);
req1[16] = MPI_COMM_SCALBL.Isend(&sendCount_Yz,1,rank_Yz,sendtag+16);
req2[16] = MPI_COMM_SCALBL.Irecv(&recvCount_yZ,1,rank_yZ,recvtag+16);
req1[17] = MPI_COMM_SCALBL.Isend(&sendCount_yZ,1,rank_yZ,sendtag+17);
req2[17] = MPI_COMM_SCALBL.Irecv(&recvCount_Yz,1,rank_Yz,recvtag+17);
/* Corners */
req1[18] = MPI_COMM_SCALBL.Isend(&sendCount_xyz,1,rank_xyz,sendtag+18);
req2[18] = MPI_COMM_SCALBL.Irecv(&recvCount_XYZ,1,rank_XYZ,recvtag+18);
req1[19] = MPI_COMM_SCALBL.Isend(&sendCount_XYz,1,rank_XYz,sendtag+19);
req2[19] = MPI_COMM_SCALBL.Irecv(&recvCount_xyZ,1,rank_xyZ,recvtag+19);
req1[20] = MPI_COMM_SCALBL.Isend(&sendCount_Xyz,1,rank_Xyz,sendtag+20);
req2[20] = MPI_COMM_SCALBL.Irecv(&recvCount_xYZ,1,rank_xYZ,recvtag+20);
req1[21] = MPI_COMM_SCALBL.Isend(&sendCount_xYz,1,rank_xYz,sendtag+21);
req2[21] = MPI_COMM_SCALBL.Irecv(&recvCount_XyZ,1,rank_XyZ,recvtag+21);
req1[22] = MPI_COMM_SCALBL.Isend(&sendCount_xyZ,1,rank_xyZ,sendtag+22);
req2[22] = MPI_COMM_SCALBL.Irecv(&recvCount_XYz,1,rank_XYz,recvtag+22);
req1[23] = MPI_COMM_SCALBL.Isend(&sendCount_XYZ,1,rank_XYZ,sendtag+23);
req2[23] = MPI_COMM_SCALBL.Irecv(&recvCount_xyz,1,rank_xyz,recvtag+23);
req1[24] = MPI_COMM_SCALBL.Isend(&sendCount_XyZ,1,rank_XyZ,sendtag+24);
req2[24] = MPI_COMM_SCALBL.Irecv(&recvCount_xYz,1,rank_xYz,recvtag+24);
req1[25] = MPI_COMM_SCALBL.Isend(&sendCount_xYZ,1,rank_xYZ,sendtag+25);
req2[25] = MPI_COMM_SCALBL.Irecv(&recvCount_Xyz,1,rank_Xyz,recvtag+25);
//...................................................................................
}
ScaLBLWideHalo_Communicator::~ScaLBLWideHalo_Communicator()
{
}
void ScaLBLWideHalo_Communicator::Recv(double *data){
//...................................................................................
MPI_Waitall(26,req1,stat1);
MPI_Waitall(26,req2,stat2);
ScaLBL_DeviceBarrier();
//...................................................................................
//...................................................................................
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_X, recvCount_X,recvbuf_X, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Y, recvCount_Y,recvbuf_Y, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xy, recvCount_xy,recvbuf_xy, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xY, recvCount_xY,recvbuf_xY, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_XY, recvCount_XY,recvbuf_XY, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_z, recvCount_z,recvbuf_z, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xz, recvCount_xz,recvbuf_xz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_yz, recvCount_yz,recvbuf_yz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Z, recvCount_Z,recvbuf_Z, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, data, Nh);
/* corners */
ScaLBL_Scalar_Unpack(dvcRecvList_xyz, recvCount_xyz,recvbuf_xyz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xYz, recvCount_xYz,recvbuf_xYz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xyZ, recvCount_xyZ,recvbuf_xyZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_xYZ, recvCount_xYZ,recvbuf_xYZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_Xyz, recvCount_Xyz,recvbuf_Xyz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_XYz, recvCount_XYz,recvbuf_XYz, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_XyZ, recvCount_XyZ,recvbuf_XyZ, data, Nh);
ScaLBL_Scalar_Unpack(dvcRecvList_XYZ, recvCount_XYZ,recvbuf_XYZ, data, Nh);
//...................................................................................
Lock=false; // unlock the communicator after communications complete
//...................................................................................
}

116
common/WideHalo.h Normal file
View File

@@ -0,0 +1,116 @@
/*
This class implements support for halo widths larger than 1
*/
#ifndef WideHalo_H
#define WideHalo_H
#include "common/ScaLBL.h"
class ScaLBLWideHalo_Communicator{
public:
//......................................................................................
ScaLBLWideHalo_Communicator(std::shared_ptr <Domain> Dm, int width);
~ScaLBLWideHalo_Communicator();
//......................................................................................
//MPI_Comm MPI_COMM_SCALBL; // MPI Communicator
Utilities::MPI MPI_COMM_SCALBL;
unsigned long int CommunicationCount,SendCount,RecvCount;
int Nx,Ny,Nz,N; // original domain structure
int Nxh,Nyh,Nzh,Nh; // with wide halo
DoubleArray Map; // map to regular halo
int first_interior,last_interior;
//......................................................................................
// Set up for D3Q19 distributions -- all 27 neighbors are needed
//......................................................................................
// Buffers to store data sent and recieved by this MPI process
double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z;
double *sendbuf_xy, *sendbuf_yz, *sendbuf_xz, *sendbuf_Xy, *sendbuf_Yz, *sendbuf_xZ;
double *sendbuf_xY, *sendbuf_yZ, *sendbuf_Xz, *sendbuf_XY, *sendbuf_YZ, *sendbuf_XZ;
double *sendbuf_xyz, *sendbuf_Xyz, *sendbuf_xYz, *sendbuf_XYz;
double *sendbuf_xyZ, *sendbuf_XyZ, *sendbuf_xYZ, *sendbuf_XYZ;
double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z;
double *recvbuf_xy, *recvbuf_yz, *recvbuf_xz, *recvbuf_Xy, *recvbuf_Yz, *recvbuf_xZ;
double *recvbuf_xY, *recvbuf_yZ, *recvbuf_Xz, *recvbuf_XY, *recvbuf_YZ, *recvbuf_XZ;
double *recvbuf_xyz, *recvbuf_Xyz, *recvbuf_xYz, *recvbuf_XYz;
double *recvbuf_xyZ, *recvbuf_XyZ, *recvbuf_xYZ, *recvbuf_XYZ;
//......................................................................................
int LastExterior();
int FirstInterior();
int LastInterior();
void Send(double *data);
void Recv(double *data);
// Debugging and unit testing functions
void PrintDebug();
private:
bool Lock; // use Lock to make sure only one call at a time to protect data in transit
// only one set of Send requests can be active at any time (per instance)
int i,j,k,n;
int iproc,jproc,kproc;
int nprocx,nprocy,nprocz;
int sendtag,recvtag;
// Give the object it's own MPI communicator
RankInfoStruct rank_info;
MPI_Group Group; // Group of processors associated with this domain
MPI_Request req1[26],req2[26];
MPI_Status stat1[26],stat2[26];
//......................................................................................
// MPI ranks for all 18 neighbors
//......................................................................................
// These variables are all private to prevent external things from modifying them!!
//......................................................................................
int rank;
int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
int rank_xy,rank_XY,rank_xY,rank_Xy;
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
int rank_xyz,rank_Xyz,rank_xYz,rank_XYz;
int rank_xyZ,rank_XyZ,rank_xYZ,rank_XYZ;
//......................................................................................
//......................................................................................
int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
int sendCount_xyz,sendCount_Xyz,sendCount_xYz,sendCount_XYz;
int sendCount_xyZ,sendCount_XyZ,sendCount_xYZ,sendCount_XYZ;
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;
int recvCount_xyz,recvCount_Xyz,recvCount_xYz,recvCount_XYz;
int recvCount_xyZ,recvCount_XyZ,recvCount_xYZ,recvCount_XYZ;
//......................................................................................
// Send buffers that reside on the compute device
int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X, *dvcSendList_Y, *dvcSendList_Z;
int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ;
int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ;
int *dvcSendList_xyz,*dvcSendList_Xyz,*dvcSendList_xYz,*dvcSendList_XYz;
int *dvcSendList_xyZ,*dvcSendList_XyZ,*dvcSendList_xYZ,*dvcSendList_XYZ;
// Recieve buffers that reside on the compute device
int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X, *dvcRecvList_Y, *dvcRecvList_Z;
int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ;
int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ;
int *dvcRecvList_xyz,*dvcRecvList_Xyz,*dvcRecvList_xYz,*dvcRecvList_XYz;
int *dvcRecvList_xyZ,*dvcRecvList_XyZ,*dvcRecvList_xYZ,*dvcRecvList_XYZ;
//......................................................................................
inline int getHaloBlock(int imin, int imax, int jmin, int jmax, int kmin, int kmax, int *dvcList){
int count = 0;
int *List;
List = new int [(imax-imin)*(jmax-jmin)*(kmax-kmin)];
for (int k=kmin; k<kmax; k++){
for (j=jmin; j<jmax; j++){
for (i=imin; i<imax; i++){
List[count++] = k*Nxh*Nyh + j*Nxh + i;
}
}
}
size_t numbytes=count*sizeof(int);
ScaLBL_AllocateZeroCopy((void **) &dvcList, numbytes); // Allocate device memory
ScaLBL_CopyToZeroCopy(dvcList,List,numbytes);
return count;
}
};
#endif

2820
cpu/FreeLee.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

48
cpu/MixedGradient.cpp Normal file
View File

@@ -0,0 +1,48 @@
/* Implement Mixed Gradient (Lee et al. JCP 2016)*/
extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz)
{
static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
{1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},
{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
int i,j,k,n,N;
int np,np2,nm; // neighbors
double v,vp,vp2,vm; // values at neighbors
double grad;
for (int idx=start; idx<finish; idx++){
n = Map[idx]; // layout in regular array
//.......Back out the 3-D indices for node n..............
k = n/(Nx*Ny);
j = (n-Nx*Ny*k)/Nx;
i = n-Nx*Ny*k-Nx*j;
v = Phi[n];
grad = 0.0;
for (int q=0; q<6; q++){
int iqx = D3Q19[q][0];
int iqy = D3Q19[q][1];
int iqz = D3Q19[q][2];
np = (k+iqz)*Nx*Ny + (j+iqy)*Nx + i + iqx;
np2 = (k+2*iqz)*Nx*Ny + (j+2*iqy)*Nx + i + 2*iqx;
nm = (k-iqz)*Nx*Ny + (j-iqy)*Nx + i - iqx;
vp = Phi[np];
vp2 = Phi[np2];
vm = Phi[nm];
grad += 0.25*(5.0*vp-vp2-3.0*v-vm);
}
for (int q=6; q<18; q++){
int iqx = D3Q19[q][0];
int iqy = D3Q19[q][1];
int iqz = D3Q19[q][2];
np = (k+iqz)*Nx*Ny + (j+iqy)*Nx + i + iqx;
np2 = (k+2*iqz)*Nx*Ny + (j+2*iqy)*Nx + i + 2*iqx;
nm = (k-iqz)*Nx*Ny + (j-iqy)*Nx + i - iqx;
vp = Phi[np];
vp2 = Phi[np2];
vm = Phi[nm];
grad += 0.125*(5.0*vp-vp2-3.0*v-vm);
}
Gradient[n] = grad;
}
}

2
models/ColorModel.cpp
View File

@@ -333,7 +333,7 @@ void ScaLBL_ColorModel::Create(){
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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................

2
models/DFHModel.cpp
View File

@@ -212,7 +212,7 @@ void ScaLBL_DFHModel::Create(){
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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
ScaLBL_Comm->Barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE

563
models/FreeLeeModel.cpp Normal file
View File

@@ -0,0 +1,563 @@
/*
color lattice boltzmann model
*/
#include "models/FreeLeeModel.h"
#include "analysis/distance.h"
#include "analysis/morphology.h"
#include "common/Communication.h"
#include "common/ReadMicroCT.h"
#include <stdlib.h>
#include <time.h>
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),rhoA(0),rhoB(0),W(0),gamma(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(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)
{
}
ScaLBL_FreeLeeModel::~ScaLBL_FreeLeeModel(){
}
void ScaLBL_FreeLeeModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
freelee_db = db->getDatabase( "FreeLee" );
analysis_db = db->getDatabase( "Analysis" );
vis_db = db->getDatabase( "Visualization" );
// set defaults
timestepMax = 100000;
tauA = tauB = 1.0;
rhoA = rhoB = 1.0;
Fx = Fy = Fz = 0.0;
gamma=1e-3;
W=5;
Restart=false;
din=dout=1.0;
flux=0.0;
// Color Model parameters
if (freelee_db->keyExists( "timestepMax" )){
timestepMax = freelee_db->getScalar<int>( "timestepMax" );
}
if (freelee_db->keyExists( "tauA" )){
tauA = freelee_db->getScalar<double>( "tauA" );
}
if (freelee_db->keyExists( "tauB" )){
tauB = freelee_db->getScalar<double>( "tauB" );
}
if (freelee_db->keyExists( "rhoA" )){
rhoA = freelee_db->getScalar<double>( "rhoA" );
}
if (freelee_db->keyExists( "rhoB" )){
rhoB = freelee_db->getScalar<double>( "rhoB" );
}
if (freelee_db->keyExists( "F" )){
Fx = freelee_db->getVector<double>( "F" )[0];
Fy = freelee_db->getVector<double>( "F" )[1];
Fz = freelee_db->getVector<double>( "F" )[2];
}
if (freelee_db->keyExists( "gamma" )){
gamma = freelee_db->getScalar<double>( "gamma" );
}
if (freelee_db->keyExists( "W" )){
W = freelee_db->getScalar<double>( "W" );
}
if (freelee_db->keyExists( "Restart" )){
Restart = freelee_db->getScalar<bool>( "Restart" );
}
if (freelee_db->keyExists( "din" )){
din = freelee_db->getScalar<double>( "din" );
}
if (freelee_db->keyExists( "dout" )){
dout = freelee_db->getScalar<double>( "dout" );
}
if (freelee_db->keyExists( "flux" )){
flux = freelee_db->getScalar<double>( "flux" );
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
outletB=1.f;
//if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
BoundaryCondition = 0;
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
}
void ScaLBL_FreeLeeModel::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
// domain parameters
Nx = Dm->Nx;
Ny = Dm->Ny;
Nz = Dm->Nz;
Lx = Dm->Lx;
Ly = Dm->Ly;
Lz = Dm->Lz;
N = Nx*Ny*Nz;
Nxh = Nx+2;
Nyh = Ny+2;
Nzh = Nz+2;
Nh = Nxh*Nyh*Nzh;
id = new signed char [N];
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
comm.barrier();
Dm->CommInit();
comm.barrier();
// Read domain parameters
rank = Dm->rank();
nprocx = Dm->nprocx();
nprocy = Dm->nprocy();
nprocz = Dm->nprocz();
}
void ScaLBL_FreeLeeModel::ReadInput(){
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
if (freelee_db->keyExists( "image_sequence" )){
auto ImageList = freelee_db->getVector<std::string>( "image_sequence");
int IMAGE_INDEX = freelee_db->getWithDefault<int>( "image_index", 0 );
std::string first_image = ImageList[IMAGE_INDEX];
Mask->Decomp(first_image);
IMAGE_INDEX++;
}
else if (domain_db->keyExists( "GridFile" )){
// Read the local domain data
auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
// Fill the halo (assuming GCW of 1)
array<int,3> size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) };
ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz };
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id.data() );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
else if (domain_db->keyExists( "Filename" )){
auto Filename = domain_db->getScalar<std::string>( "Filename" );
Mask->Decomp(Filename);
}
else{
Mask->ReadIDs();
}
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz);
// Solve for the position of the solid phase
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
// Initialize the solid phase
signed char label = Mask->id[n];
if (label > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
}
}
}
SignDist.resize(Nx,Ny,Nz);
// Initialize the signed distance function
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
// Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(SignDist,id_solid,*Mask);
if (rank == 0) cout << "Domain set." << endl;
}
void ScaLBL_FreeLeeModel::Create(){
/*
* This function creates the variables needed to run a LBM
*/
//.........................................................
// 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,2);
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 **) &fq, 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] = k*Nx*Ny+j*Nx+i;
}
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
auto n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n", n);
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
auto n = TmpMap[idx];
if ( n > Nx*Ny*Nz ){
printf("Bad value! idx=%i \n",n);
TmpMap[idx] = Nx*Ny*Nz-1;
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
delete [] TmpMap;
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
// initialize phi based on PhaseLabel (include solid component labels)
}
/********************************************************
* AssignComponentLabels *
********************************************************/
void ScaLBL_FreeLeeModel::Initialize(){
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
/*
* This function initializes model
*/
if (Restart == true){
if (rank==0){
printf("Reading restart file! \n");
}
// Read in the restart file to CPU buffers
int *TmpMap;
TmpMap = new int[Np];
double *cPhi, *cDist, *cDen;
cPhi = new double[N];
cDen = new double[2*Np];
cDist = new double[19*Np];
ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int));
ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double));
ifstream File(LocalRestartFile,ios::binary);
int idx;
double value,va,vb;
for (int n=0; n<Np; n++){
File.read((char*) &va, sizeof(va));
File.read((char*) &vb, sizeof(vb));
cDen[n] = va;
cDen[Np+n] = vb;
}
for (int n=0; n<Np; n++){
// Read the distributions
for (int q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cDist[q*Np+n] = value;
}
}
File.close();
for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
va = cDen[n];
vb = cDen[Np + n];
value = (va-vb)/(va+vb);
idx = TmpMap[n];
if (!(idx < 0) && idx<N)
cPhi[idx] = value;
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
va = cDen[n];
vb = cDen[Np + n];
value = (va-vb)/(va+vb);
idx = TmpMap[n];
if (!(idx < 0) && idx<N)
cPhi[idx] = value;
}
// Copy the restart data to the GPU
ScaLBL_CopyToDevice(Den,cDen,2*Np*sizeof(double));
ScaLBL_CopyToDevice(fq,cDist,19*Np*sizeof(double));
ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
ScaLBL_DeviceBarrier();
comm.barrier();
}
if (rank==0) printf ("Initializing phase field \n");
//ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_PhaseField_Init(dvcMap, Phi, Den, hq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
// establish reservoirs for external bC
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4 ){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
//ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
}
void ScaLBL_FreeLeeModel::Run(){
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);
}
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier();
comm.barrier();
starttime = MPI_Wtime();
//.........................................
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
while (timestep < timestepMax ) {
//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 hq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(hq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(hq,Bq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, 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, hq, 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_DeviceBarrier();
// 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, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(hq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(hq,Bq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, 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, hq, 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_DeviceBarrier();
// 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, hq, 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);
//************************************************************************
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
// ************************************************************************
}
void ScaLBL_FreeLeeModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);
//ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
ScaLBL_CopyToHost(PhaseField.data(), Phi, sizeof(double)*N);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank);
AFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,AFILE);
fclose(AFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
FILE *BFILE;
sprintf(LocalRankFilename,"B.%05i.raw",rank);
BFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,BFILE);
fclose(BFILE);
ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
FILE *PFILE;
sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
PFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,PFILE);
fclose(PFILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
VELX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELX_FILE);
fclose(VELX_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
FILE *VELY_FILE;
sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
VELY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELY_FILE);
fclose(VELY_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
FILE *VELZ_FILE;
sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
VELZ_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE);
/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
FILE *CGX_FILE;
sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
CGX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGX_FILE);
fclose(CGX_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField);
FILE *CGY_FILE;
sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank);
CGY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGY_FILE);
fclose(CGY_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField);
FILE *CGZ_FILE;
sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank);
CGZ_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGZ_FILE);
fclose(CGZ_FILE);
*/
}

86
models/FreeLeeModel.h Normal file
View File

@@ -0,0 +1,86 @@
/*
Implementation of Lee et al JCP 2016 lattice boltzmann model
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/Communication.h"
#include "common/MPI.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
#include "common/ScaLBL.h"
#include "common/WideHalo.h"
class ScaLBL_FreeLeeModel{
public:
ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM);
~ScaLBL_FreeLeeModel();
// functions in they should be run
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void Run();
void WriteDebug();
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
double tauA,tauB,rhoA,rhoB;
double W,gamma;
double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB;
int Nx,Ny,Nz,N,Np;
int Nxh,Nyh,Nzh,Nh; // extra halo width
int rank,nprocx,nprocy,nprocz,nprocs;
double Lx,Ly,Lz;
std::shared_ptr<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
std::shared_ptr<ScaLBLWideHalo_Communicator> ScaLBL_Comm_WideHalo;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
std::shared_ptr<Database> freelee_db;
std::shared_ptr<Database> analysis_db;
std::shared_ptr<Database> vis_db;
IntArray Map;
signed char *id;
int *NeighborList;
int *dvcMap;
double *fq, *hq;
double *mu_phi, *Den, *Phi;
double *ColorGrad;
double *Velocity;
double *Pressure;
DoubleArray SignDist;
private:
Utilities::MPI comm;
int dist_mem_size;
int neighborSize;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
};

2
models/GreyscaleColorModel.cpp
View File

@@ -600,7 +600,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................

2
models/GreyscaleModel.cpp
View File

@@ -366,7 +366,7 @@ void ScaLBL_GreyscaleModel::Create(){
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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................

3
models/IonModel.cpp
View File

@@ -650,8 +650,9 @@ void ScaLBL_IonModel::Create(){
if (rank==0) printf ("LB Ion Solver: Set up memory efficient layout \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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................

3
models/MRTModel.cpp
View File

@@ -172,8 +172,9 @@ void ScaLBL_MRTModel::Create(){
if (rank==0) printf ("Set up memory efficient layout \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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................

3
models/PoissonSolver.cpp
View File

@@ -275,8 +275,9 @@ void ScaLBL_Poisson::Create(){
if (rank==0) printf ("LB-Poisson Solver: Set up memory efficient layout \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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................

3
models/StokesModel.cpp
View File

@@ -278,8 +278,9 @@ void ScaLBL_StokesModel::Create(){
if (rank==0) printf ("LB Single-Fluid Solver: Set up memory efficient layout \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);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................

13
tests/CMakeLists.txt
View File

@@ -35,8 +35,13 @@ ADD_LBPM_EXECUTABLE( GenerateSphereTest )
#ADD_LBPM_EXECUTABLE( BlobIdentifyParallel )
#ADD_LBPM_EXECUTABLE( convertIO )
#ADD_LBPM_EXECUTABLE( DataAggregator )
#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )
#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )(
ADD_LBPM_EXECUTABLE( lbpm_minkowski_scalar )
ADD_LBPM_EXECUTABLE( TestPoissonSolver )
ADD_LBPM_EXECUTABLE( TestIonModel )
ADD_LBPM_EXECUTABLE( TestNernstPlanck )
ADD_LBPM_EXECUTABLE( TestPNP_Stokes )
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY )
@@ -49,13 +54,9 @@ ADD_LBPM_TEST( TestTorusEvolve )
ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestPoissonSolver )
ADD_LBPM_TEST( TestIonModel )
ADD_LBPM_TEST( TestNernstPlanck )
ADD_LBPM_TEST( TestPNP_Stokes )
#ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad )
#ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestWideHalo )
ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db)
#ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db)

2
tests/TestBubbleDFH.cpp
View File

@@ -247,7 +247,7 @@ int main(int argc, char **argv)
if (rank==0) printf ("Set up memory efficient layout Npad=%i \n",Npad);
IntArray Map(Nx,Ny,Nz);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1);
comm.barrier();
//...........................................................................

107
tests/TestColorGrad.cpp
View File

@@ -3,11 +3,16 @@
// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media
// James E. McCLure
//*************************************************************************
#include <stdio.h>
#include <stdio.h> // Initialize MPI
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
int check;
#include <iostream>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
using namespace std;
@@ -47,105 +52,13 @@ int main(int argc, char **argv)
int Nx,Ny,Nz;
int i,j,k,n;
int dim = 3;
Nx = Ny = Nz = 32;
Lx = Ly = Lz = 1.0;
//if (rank == 0) printf("dim=%d\n",dim);
int timestep = 0;
int timesteps = 100;
int centralNode = 2;
double tauA = 1.0;
double tauB = 1.0;
double rhoA = 1.0;
double rhoB = 1.0;
double alpha = 0.005;
double beta = 0.95;
double tau = 1.0;
double mu=(tau-0.5)/3.0;
double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
Fx = Fy = 0.f;
Fz = 0.f;
if (rank==0){
//.......................................................................
// Reading the domain information file
//.......................................................................
ifstream domain("Domain.in");
if (domain.good()){
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> Nx;
domain >> Ny;
domain >> Nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
}
else if (nprocs==1){
nprocx=nprocy=nprocz=1;
Nx=Ny=Nz=3;
nspheres=0;
Lx=Ly=Lz=1;
}
else if (nprocs==2){
nprocx=2; nprocy=1;
nprocz=1;
Nx=Ny=Nz=dim;
Nx = dim; Ny = dim; Nz = dim;
nspheres=0;
Lx=Ly=Lz=1;
}
else if (nprocs==4){
nprocx=nprocy=2;
nprocz=1;
Nx=Ny=Nz=dim;
nspheres=0;
Lx=Ly=Lz=1;
}
else if (nprocs==8){
nprocx=nprocy=nprocz=2;
Nx=Ny=Nz=dim;
nspheres=0;
Lx=Ly=Lz=1;
}
//.......................................................................
}
// **************************************************************
// Broadcast simulation parameters from rank 0 to all other procs
MPI_Barrier(comm);
//.................................................
MPI_Bcast(&Nx,1,MPI_INT,0,comm);
MPI_Bcast(&Ny,1,MPI_INT,0,comm);
MPI_Bcast(&Nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
MPI_Barrier(comm);
// **************************************************************
// **************************************************************
if (nprocs != nprocx*nprocy*nprocz){
printf("nprocx = %i \n",nprocx);
printf("nprocy = %i \n",nprocy);
printf("nprocz = %i \n",nprocz);
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
if (rank==0){
printf("********************************************************\n");
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("********************************************************\n");
}
MPI_Barrier(comm);
double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz;
int BoundaryCondition=0;
@@ -190,7 +103,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np,1);
MPI_Barrier(comm);
//......................device distributions.................................

2
tests/TestColorGradDFH.cpp
View File

@@ -100,7 +100,7 @@ int main(int argc, char **argv)
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
//......................device distributions.................................

2
tests/TestColorMassBounceback.cpp
View File

@@ -168,7 +168,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
Npad=Np+32;
neighborList= new int[18*Npad];
Np=ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
Np=ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np,1);
MPI_Barrier(comm);
//......................device distributions.................................

3
tests/TestCommD3Q19.cpp
View File

@@ -283,8 +283,9 @@ int main(int argc, char **argv)
auto neighborList= new int[18*Npad];
IntArray Map(Nx,Ny,Nz);
Map.fill(-2);
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
int neighborSize=18*Np*sizeof(int);
//......................device distributions.................................
dist_mem_size = Np*sizeof(double);

3
tests/TestFluxBC.cpp
View File

@@ -88,8 +88,9 @@ int main (int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
if (rank==0) printf ("Allocating distributions \n");

3
tests/TestForceMoments.cpp
View File

@@ -179,8 +179,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
if (rank == 0) PrintNeighborList(neighborList,Np, rank);

7
tests/TestMRT.cpp
View File

@@ -704,7 +704,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np,1);
MPI_Barrier(comm);
//......................device distributions.................................
@@ -804,13 +804,8 @@ int main(int argc, char **argv)
}
// ****************************************************
<<<<<<< HEAD
comm.barrier();
Utilities::shutdown();
=======
MPI_Barrier(comm);
MPI_Finalize();
>>>>>>> electrokinetic
// ****************************************************
return check;

2
tests/TestMap.cpp
View File

@@ -88,7 +88,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
// Check the neighborlist

2
tests/TestPressVel.cpp
View File

@@ -127,7 +127,7 @@ int main(int argc, char **argv)
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
//......................device distributions.................................

204
tests/TestWideHalo.cpp Normal file
View File

@@ -0,0 +1,204 @@
//*************************************************************************
// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media
// James E. McCLure
//*************************************************************************
#include <stdio.h>
#include <iostream>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/WideHalo.h"
#include "common/MPI.h"
using namespace std;
//***************************************************************************************
int main(int argc, char **argv)
{
//*****************************************
// ***** MPI STUFF ****************
//*****************************************
// Initialize MPI
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
int check=0;
{
if (rank == 0){
printf("********************************************************\n");
printf("Running Color Model: TestColor \n");
printf("********************************************************\n");
}
// Domain variables
int nprocx, nprocy, nprocz;
double Lx,Ly,Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int dim = 16;
Lx = Ly = Lz = 1.0;
int BoundaryCondition=0;
//.......................................................................
// Reading the domain information file
//.......................................................................
nprocx=nprocy=nprocz=1;
if (nprocs==1){
nprocx=nprocy=nprocz=1;
Nx=Ny=Nz=dim;
Lx=Ly=Lz=1;
}
else if (nprocs==2){
nprocx=2; nprocy=1;
nprocz=1;
Nx=Ny=Nz=dim;
Nx = dim; Ny = dim; Nz = dim;
Lx=Ly=Lz=1;
}
else if (nprocs==4){
nprocx=nprocy=2;
nprocz=1;
Nx=Ny=Nz=dim;
Lx=Ly=Lz=1;
}
else if (nprocs==8){
nprocx=nprocy=nprocz=2;
Nx=Ny=Nz=dim;
Lx=Ly=Lz=1;
}
//.......................................................................
// **************************************************************
if (nprocs != nprocx*nprocy*nprocz){
printf("nprocx = %i \n",nprocx);
printf("nprocy = %i \n",nprocy);
printf("nprocz = %i \n",nprocz);
INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
}
if (rank==0){
printf("********************************************************\n");
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("********************************************************\n");
}
comm.barrier();
std::shared_ptr<Domain> Dm = std::shared_ptr<Domain>(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition));
Nx += 2;
Ny += 2;
Nz += 2;
int N = Nx*Ny*Nz;
int Np=0; // number of local pore nodes
double *PhaseLabel;
PhaseLabel = new double[N];
//.......................................................................
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Dm->id[n]=1;
Np++;
// Initialize gradient ColorGrad = (1,2,3)
double value=double(3*k+2*j+i);
PhaseLabel[n]= value;
}
}
}
Dm->CommInit();
comm.barrier();
if (rank == 0) cout << "Domain set." << endl;
if (rank==0) printf ("Create ScaLBL_Communicator \n");
//Create a second communicator based on the regular data layout
ScaLBL_Communicator ScaLBL_Comm_Regular(Dm);
ScaLBL_Communicator ScaLBL_Comm(Dm);
ScaLBLWideHalo_Communicator WideHalo(Dm,2);
// LBM variables
if (rank==0) printf ("Set up the neighborlist \n");
int neighborSize=18*Np*sizeof(int);
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,2);
comm.barrier();
int *NeighborList;
int *dvcMap;
double *Phi;
double *ColorGrad;
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
int *WideMap;
WideMap=new int[Np];
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
int nw = WideHalo.Map(i,j,k);
int idx = Map(i,j,k);
WideMap[idx] = nw;
}
}
}
ScaLBL_CopyToDevice(dvcMap, WideMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
// initialize phi based on PhaseLabel (include solid component labels)
ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double));
//...........................................................................
int Nxh = Nx+2;
int Nyh = Ny+2;
int Nzh = Nz+2;
ScaLBL_D3Q19_MixedGradient(dvcMap, Phi, ColorGrad, 0, Np, Np, Nxh, Nyh, Nzh);
double *COLORGRAD;
COLORGRAD= new double [3*Np];
int SIZE=3*Np*sizeof(double);
ScaLBL_CopyToHost(&COLORGRAD[0],&ColorGrad[0],SIZE);
double CX,CY,CZ;
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
CX=COLORGRAD[idx];
CY=COLORGRAD[Np+idx];
CZ=COLORGRAD[2*Np+idx];
double error=sqrt((CX-1.0)*(CX-1.0)+(CY-2.0)*(CY-2.0)+ (CZ-3.0)*(CZ-3.0));
if (error > 1e-8){
check++;
printf("i,j,k=%i,%i,%i: Color gradient=%f,%f,%f \n",i,j,k,CX,CY,CZ);
}
}
}
}
}
}
// ****************************************************
comm.barrier();
Utilities::shutdown();
// ****************************************************
return check;
}