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

adding mass conservation test

Browse Source
This commit is contained in:
James McClure 2021-03-11 13:55:47 -05:00
parent 9c04320276
commit 3f659a0cc2
3 changed files with 114 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
cpu/FreeLee.cpp
View File

@ -778,25 +778,25 @@ extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel(int *neighborList, int *Map, dou
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
hq[n] = h0 - (h0 - 0.33333333333333333*phi)/tauM;
// q = 1
hq[nr2] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
hq[nr2] = h1 - (h1 - 0.1111111111111111*(phi*(1.0 + 3.0*ux) - (M*nx*(1 - 4*phi*phi))/W)/tauM);
// q = 2
hq[nr1] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
hq[nr1] = h2 - (h2 - 0.1111111111111111*(phi*(1.0 - 3.0*ux) + (M*nx*(1 - 4*phi*phi))/W)/tauM);
// q = 3
hq[nr4] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
hq[nr4] = h3 - (h3 - 0.1111111111111111*(phi*(1.0 + 3.0*uy) - (M*ny*(1 - 4*phi*phi))/W)/tauM);
// q = 4
hq[nr3] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
hq[nr3] = h4 - (h4 - 0.1111111111111111*(phi*(1.0 - 3.0*uy) + (M*ny*(1 - 4*phi*phi))/W)/tauM);
// q = 5
hq[nr6] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
hq[nr6] = h5 - (h5 - 0.1111111111111111*(phi*(1.0 + 3.0*uz) - (M*nz*(1 - 4*phi*phi))/W)/tauM);
// q = 6
hq[nr5] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
hq[nr5] = h6 - (h6 - 0.1111111111111111*(phi*(1.0 - 3.0*uz) + (M*nz*(1 - 4*phi*phi))/W)/tauM);
//........................................................................
//Update velocity on device
@ -1369,25 +1369,25 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
//-------------------------------- BGK collison for phase field ---------------------------------//
// q = 0
hq[n] = h0 - (h0 - 0.3333333333333333*phi)/tauM;
hq[n] = h0 - (h0 - 0.33333333333333333*phi)/tauM;
// q = 1
hq[1*Np+n] = h1 - (h1 - phi*(0.1111111111111111 + 0.5*ux) - (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
hq[1*Np+n] = h1 - (h1 - 0.1111111111111111*(phi*(1.0 + 3.0*ux) - (M*nx*(1 - 4*phi*phi))/W)/tauM);
// q = 2
hq[2*Np+n] = h2 - (h2 - phi*(0.1111111111111111 - 0.5*ux) + (0.5*M*nx*(1 - 4*phi*phi))/W)/tauM;
hq[2*Np+n] = h2 - (h2 - 0.1111111111111111*(phi*(1.0 - 3.0*ux) + (M*nx*(1 - 4*phi*phi))/W)/tauM);
// q = 3
hq[3*Np+n] = h3 - (h3 - phi*(0.1111111111111111 + 0.5*uy) - (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
hq[3*Np+n] = h3 - (h3 - 0.1111111111111111*(phi*(1.0 + 3.0*uy) - (M*ny*(1 - 4*phi*phi))/W)/tauM);
// q = 4
hq[4*Np+n] = h4 - (h4 - phi*(0.1111111111111111 - 0.5*uy) + (0.5*M*ny*(1 - 4*phi*phi))/W)/tauM;
hq[4*Np+n] = h4 - (h4 - 0.1111111111111111*(phi*(1.0 - 3.0*uy) + (M*ny*(1 - 4*phi*phi))/W)/tauM);
// q = 5
hq[5*Np+n] = h5 - (h5 - phi*(0.1111111111111111 + 0.5*uz) - (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
hq[5*Np+n] = h5 - (h5 - 0.1111111111111111*(phi*(1.0 + 3.0*uz) - (M*nz*(1 - 4*phi*phi))/W)/tauM);
// q = 6
hq[6*Np+n] = h6 - (h6 - phi*(0.1111111111111111 - 0.5*uz) + (0.5*M*nz*(1 - 4*phi*phi))/W)/tauM;
hq[6*Np+n] = h6 - (h6 - 0.1111111111111111*(phi*(1.0 - 3.0*uz) + (M*nz*(1 - 4*phi*phi))/W)/tauM);
//........................................................................
//Update velocity on device

1
tests/CMakeLists.txt
View File

@ -62,6 +62,7 @@ ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestWideHalo )
ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db)
ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
#ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db)

99
tests/testGlobalMassFreeLee.cpp Normal file
View File

@ -0,0 +1,99 @@
#include <exception>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include "common/Utilities.h"
#include "models/FreeLeeModel.h"
//*******************************************************************
// Implementation of Free-Energy Two-Phase LBM (Lee model)
//*******************************************************************
int main( int argc, char **argv )
{
// Initialize
Utilities::startup( argc, argv );
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
if (rank == 0){
printf("********************************************************\n");
printf("Running Free Energy Lee LBM \n");
printf("********************************************************\n");
}
// Initialize compute device
int device=ScaLBL_SetDevice(rank);
NULL_USE( device );
ScaLBL_DeviceBarrier();
comm.barrier();
PROFILE_ENABLE(1);
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_FreeLeeModel LeeModel( rank,nprocs,comm );
LeeModel.ReadParams( filename );
LeeModel.SetDomain();
LeeModel.ReadInput();
LeeModel.Create_TwoFluid();
LeeModel.Initialize_TwoFluid();
/* Copy the initial density to test that global mass is conserved */
int Nx = LeeModel.Dm->Nx;
int Ny = LeeModel.Dm->Ny;
int Nz = LeeModel.Dm->Nz;
DoubleArray DensityInit(Nx,Ny,Nz);
LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityInit);
LeeModel.Run_TwoFluid();
DoubleArray DensityFinal(Nx,Ny,Nz);
LeeModel.ScaLBL_Comm->RegularLayout(LeeModel.Map,LeeModel.Den,DensityFinal);
DoubleArray DensityChange(Nx,Ny,Nz);
double totalChange=0.0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
double change = DensityFinal(i,j,k)-DensityInit(i,j,k);
DensityChange(i,j,k) = change;
totalChange += change;
}
}
}
printf("Density change, %f\n", totalChange);
FILE *OUTFILE;
char LocalRankFilename[40];
sprintf(LocalRankFilename,"DensityChange.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(DensityChange.data(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
LeeModel.WriteDebug_TwoFluid();
PROFILE_STOP("Main");
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
return 0;
}