add a test routine for mixed gradient

This commit is contained in:
Rex Zhe Li
2021-02-25 20:30:48 -05:00
parent a0b42380a4
commit 99f1d9b727
6 changed files with 402 additions and 0 deletions

View File

@@ -207,6 +207,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,
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 start, int finish, Np);
// BOUNDARY CONDITION ROUTINES

View File

@@ -1916,3 +1916,167 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, d
Pressure[n] = p;
}
}
extern "C" void ScaLBL_D3Q9_MGTest(int *Map, double *Phi,double *ColorGrad, int start, int finish, Np){
int n,nn,nn2x,ijk;
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,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;
}
}

View File

@@ -1059,3 +1059,151 @@ void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){
fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE);
}
void ScaLBL_FreeLeeModel::Create_DummyPhase_MGTest(){
// Initialize communication structures in averaging domain
for (int i=0; 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 **) &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, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_D3Q9_MGTest(dvcMap,Phi,ColorGrad, 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",cgx);
printf("Sum of all y-component of the mixed gradient = %.2g",cgy);
printf("Sum of all z-component of the mixed gradient = %.2g",cgz);
}
delete [] ColorGrad_host;
}

View File

@@ -34,6 +34,9 @@ public:
void Initialize_SingleFluid();
void Run_SingleFluid();
void WriteDebug_SingleFluid();
// test utilities
void Create_DummyPhase_MGTest();
void MGTest();
bool Restart,pBC;
int timestep,timestepMax;

View File

@@ -43,6 +43,7 @@ ADD_LBPM_EXECUTABLE( TestPoissonSolver )
ADD_LBPM_EXECUTABLE( TestIonModel )
ADD_LBPM_EXECUTABLE( TestNernstPlanck )
ADD_LBPM_EXECUTABLE( TestPNP_Stokes )
ADD_LBPM_EXECUTABLE( TestMixedGrad )

85
tests/TestMixedGrad.cpp Normal file
View File

@@ -0,0 +1,85 @@
#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_DummyPhaseField(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) cout << "Setting up dummy phase field..." << endl;
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 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 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_DummyPhaseField(LeeModel);
LeeModel.Create_DummyPhase_MGTest();
LeeModel.MGTest()
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "TestMixedGrad" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
PROFILE_SAVE( file, level );
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
return 0;
}