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

debugging strange mass conservation issue

Browse Source
This commit is contained in:
James McClure
2020-04-10 15:03:15 -04:00
parent cfa40bdcba
commit bdf8539f40
5 changed files with 143 additions and 134 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
analysis/SubPhase.cpp
View File

@@ -280,7 +280,7 @@ void SubPhase::Basic(){
dir_y = 0.0;
dir_z = 1.0;
}
if (Dm->BoundaryCondition > 0 ){
if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){
// compute the pressure drop
double pressure_drop = (Pressure(Nx*Ny + Nx + 1) - 1.0) / 3.0;
double length = ((Nz-2)*Dm->nprocz());
@@ -376,8 +376,8 @@ void SubPhase::Full(){
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
if (Dm->BoundaryCondition > 0 && Dm->BoundaryCondition != 5 && Dm->kproc() == 0) kmin=4;
if (Dm->BoundaryCondition > 0 && Dm->BoundaryCondition != 5 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
imin=jmin=1;
// If inlet layers exist use these as default

14
common/ScaLBL.cpp
View File

@@ -1286,7 +1286,7 @@ void ScaLBL_Communicator::BiRecvD3Q7AA(double *Aq, double *Bq){
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,Bq,N);
}
if (BoundaryCondition == 5){
/* if (BoundaryCondition == 5){
if (kproc == 0){
ScaLBL_D3Q7_Reflection_BC_z(dvcSendList_z, Aq, sendCount_z, N);
ScaLBL_D3Q7_Reflection_BC_z(dvcSendList_z, Bq, sendCount_z, N);
@@ -1296,6 +1296,7 @@ void ScaLBL_Communicator::BiRecvD3Q7AA(double *Aq, double *Bq){
ScaLBL_D3Q7_Reflection_BC_Z(dvcSendList_Z, Bq, sendCount_Z, N);
}
}
*/
//...................................................................................
Lock=false; // unlock the communicator after communications complete
//...................................................................................
@@ -1543,10 +1544,15 @@ void ScaLBL_Communicator::RecvHalo(double *data){
ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, data, N);
ScaLBL_Scalar_Unpack(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, data, N);
}
//...................................................................................
Lock=false; // unlock the communicator after communications complete
//...................................................................................
if (BoundaryCondition == 5 && kproc == 0){
ScaLBL_CopySlice_z(data,Nx,Ny,Nz,1,0);
}
if (BoundaryCondition == 5 && kproc == nprocz-1){
ScaLBL_CopySlice_z(data,Nx,Ny,Nz,Nz-2,Nz-1);
}
}
void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, DoubleArray &regdata){
@@ -1583,7 +1589,7 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB){
if (kproc == 0) {
if (BoundaryCondition == 5){
ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,1,0);
//ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,1,0);
}
else {
// Set the phase indicator field and density on the z inlet
@@ -1596,7 +1602,7 @@ void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double
void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB){
if (kproc == nprocz-1){
if (BoundaryCondition == 5){
ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,Nz-2,Nz-1);
//ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,Nz-2,Nz-1);
}
else {
// Set the phase indicator field and density on the Z outlet

49
models/ColorModel.cpp
View File

@@ -470,7 +470,8 @@ void ScaLBL_ColorModel::Initialize(){
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
// 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);
@@ -743,7 +744,7 @@ void ScaLBL_ColorModel::Run(){
//************************************************************************
PROFILE_STOP("Update");
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){
printf("%i %f \n",timestep,din);
}
// Run the analysis
@@ -1159,7 +1160,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
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);
@@ -1447,7 +1448,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
// 7. Re-initialize phase field and density
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
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);
@@ -1516,25 +1517,25 @@ void ScaLBL_ColorModel::WriteDebug(){
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);
/* 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);
*/
}

206
tests/TestMassConservationD3Q7.cpp
View File

@@ -69,10 +69,11 @@ int main(int argc, char **argv)
// Initialize MPI
int rank,nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
// parallel domain size (# of sub-domains)
int CleanCheck = 0;
if (rank == 0){
printf("********************************************************\n");
@@ -84,68 +85,68 @@ int main(int argc, char **argv)
}
}
{
auto filename = argv[1];
ScaLBL_ColorModel CM(rank,nprocs,comm);
CM.ReadParams(filename);
CM.SetDomain();
int i,j,k,n;
int Nx,Ny,Nz,N,Np;
Nx = CM.Nx;
Ny = CM.Ny;
Nz = CM.Nz;
N = Nx*Ny*Nz;
auto filename = argv[1];
ScaLBL_ColorModel CM(rank,nprocs,comm);
CM.ReadParams(filename);
CM.SetDomain();
int i,j,k,n;
int Nx,Ny,Nz,N,Np;
Nx = CM.Nx;
Ny = CM.Ny;
Nz = CM.Nz;
N = Nx*Ny*Nz;
//CM.ReadInput();
double radius=0.4*double(Nx);
InitializeBubble(CM,radius);
CM.Create(); // creating the model will create data structure to match the pore structure and allocate variables
CM.Initialize(); // initializing the model will set initial conditions for variables
//CM.Run();
//CM.WriteDebug();
//CM.ReadInput();
double radius=0.4*double(Nx);
InitializeBubble(CM,radius);
CM.Create(); // creating the model will create data structure to match the pore structure and allocate variables
CM.Initialize(); // initializing the model will set initial conditions for variables
//CM.Run();
//CM.WriteDebug();
CM.timestepMax = 10;
CM.Run();
CM.timestepMax = 10;
CM.Run();
Np = CM.Np;
double *DenOriginal, *DenFinal;
DenOriginal = new double [2*Np];
DenFinal = new double [2*Np];
Np = CM.Np;
double *DenOriginal, *DenFinal;
DenOriginal = new double [2*Np];
DenFinal = new double [2*Np];
// Run the odd timestep
ScaLBL_CopyToHost(DenOriginal,CM.Den,2*Np*sizeof(double));
/*
// Run the odd timestep
ScaLBL_CopyToHost(DenOriginal,CM.Den,2*Np*sizeof(double));
/*
CM.ScaLBL_Comm->BiSendD3Q7AA(CM.Aq,CM.Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(CM.NeighborList, CM.dvcMap, CM.Aq, CM.Bq, CM.Den, CM.Phi, CM.ScaLBL_Comm->FirstInterior(), CM.ScaLBL_Comm->LastInterior(), CM.Np);
CM.ScaLBL_Comm->BiRecvD3Q7AA(CM.Aq,CM.Bq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_D3Q7_AAodd_PhaseField(CM.NeighborList, CM.dvcMap, CM.Aq, CM.Bq, CM.Den, CM.Phi, 0, CM.ScaLBL_Comm->LastExterior(), CM.Np);
*/
*/
CM.timestepMax = 2;
CM.Run();
int D3Q7[7][3]={{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
// Compare and make sure mass is conserved at every lattice site
auto Error = new double[N];
auto A_q = new double[7*Np];
//auto B_q = new double[7*Np];
bool CleanCheck = true;
double original,final, sum_q;
double total_mass_A_0 = 0.0;
double total_mass_B_0= 0.0;
double total_mass_A_1 = 0.0;
double total_mass_B_1= 0.0;
int count_negative_A = 0;
int count_negative_B = 0;
ScaLBL_CopyToHost(DenFinal,CM.Den,2*Np*sizeof(double));
ScaLBL_CopyToHost(A_q,CM.Aq,7*Np*sizeof(double));
for (i=0; i<N; i++) Error[i]=0.0;
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;
int idx = CM.Map(i,j,k);
if (idx < Np && idx>-1){
//printf("idx=%i\n",idx);
CM.timestepMax = 2;
CM.timestep = 0;
CM.Run();
int D3Q7[7][3]={{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
// Compare and make sure mass is conserved at every lattice site
auto Error = new double[N];
auto A_q = new double[7*Np];
//auto B_q = new double[7*Np];
double original,final, sum_q;
double total_mass_A_0 = 0.0;
double total_mass_B_0= 0.0;
double total_mass_A_1 = 0.0;
double total_mass_B_1= 0.0;
int count_negative_A = 0;
int count_negative_B = 0;
ScaLBL_CopyToHost(DenFinal,CM.Den,2*Np*sizeof(double));
ScaLBL_CopyToHost(A_q,CM.Aq,7*Np*sizeof(double));
for (i=0; i<N; i++) Error[i]=0.0;
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;
int idx = CM.Map(i,j,k);
if (idx < Np && idx>-1){
//printf("idx=%i\n",idx);
final = DenFinal[idx];
if (final < 0.0) count_negative_A++;
original = DenOriginal[idx];
@@ -153,60 +154,61 @@ int main(int argc, char **argv)
total_mass_A_1 += final;
sum_q = A_q[idx];
for (int q=1; q<7; q++){
int Cqx = D3Q7[q][0];
int Cqy = D3Q7[q][1];
int Cqz = D3Q7[q][2];
int iq = CM.Map(i-Cqx,j-Cqy,k-Cqz);
if (iq < Np && iq > -1){
sum_q += A_q[q*Np+iq];
}
else if (q%2==0){
sum_q += A_q[(q-1)*Np+idx];
}
else{
sum_q += A_q[(q+1)*Np+idx];
}
int Cqx = D3Q7[q][0];
int Cqy = D3Q7[q][1];
int Cqz = D3Q7[q][2];
int iq = CM.Map(i-Cqx,j-Cqy,k-Cqz);
if (iq < Np && iq > -1){
sum_q += A_q[q*Np+iq];
}
else if (q%2==0){
sum_q += A_q[(q-1)*Np+idx];
}
else{
sum_q += A_q[(q+1)*Np+idx];
}
}
Error[n] = sum_q - original;
/*if (fabs(DenFinal[idx] - DenOriginal[idx]) > 1e-15){
//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");
printf("Mass not conserved: WP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
CleanCheck=false;
Error[n] += final-original;
}*/
if (fabs(DenFinal[idx] - DenOriginal[idx]) > 1e-15){
//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");
//printf("Mass not conserved: WP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
CleanCheck=false;
Error[n] += final-original;
}
final = DenFinal[Np+idx];
if (final < 0.0) count_negative_B++;
original = DenOriginal[Np+idx];
total_mass_B_0 += original;
total_mass_B_1 += final;
/*if (fabs(DenFinal[Np+idx] - DenOriginal[Np+idx]) > 1e-15){
//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");
printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
CleanCheck=false;
Error[n] += final-original;
}*/
if (fabs(DenFinal[Np+idx] - DenOriginal[Np+idx]) > 1e-15){
//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");
//printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
CleanCheck=false;
Error[n] += final-original;
}
}
}
}
}
}
printf("Negative density values for A = %i \n",count_negative_A);
printf("Negative density values for B = %i \n",count_negative_B);
printf("Global mass difference A = %.5g\n",total_mass_A_1-total_mass_A_0);
printf("Global mass difference B = %.5g\n",total_mass_B_1-total_mass_B_0);
printf("Negative density values for A = %i \n",count_negative_A);
printf("Negative density values for B = %i \n",count_negative_B);
printf("Global mass difference A = %.5g\n",total_mass_A_1-total_mass_A_0);
printf("Global mass difference B = %.5g\n",total_mass_B_1-total_mass_B_0);
if (count_negative_A > 0 ||count_negative_B > 0) CleanCheck=1;
if (fabs(total_mass_A_1-total_mass_A_0) > 1.0e-15||fabs(total_mass_B_1-total_mass_B_0) > 1.0e-15 ) CleanCheck=2;
if (count_negative_A > 0 ||count_negative_B > 0) CleanCheck=1;
if (fabs(total_mass_A_1-total_mass_A_0) > 1.0e-13||fabs(total_mass_B_1-total_mass_B_0) > 1.0e-13 ) CleanCheck=2;
/*
FILE *OUTFILE;
OUTFILE = fopen("error.raw","wb");
fwrite(Error,8,N,OUTFILE);
fclose(OUTFILE);
FILE *OUTFILE;
OUTFILE = fopen("error.raw","wb");
fwrite(Error,8,N,OUTFILE);
fclose(OUTFILE);
/*
if (rank==0) printf("Checking that the correct velocity is retained \n");
// Swap convention is observed -- velocity is negative
@@ -256,15 +258,15 @@ int main(int argc, char **argv)
}
}
}
*/
if (CleanCheck){
if (rank==0) printf("Test passed: mass conservation for D3Q7 \n");
}
else {
if (rank==0) printf("Test failed!: mass conservation for D3Q7 \n");
*/
if (CleanCheck == 0){
if (rank==0) printf("Test passed: mass conservation for D3Q7 \n");
}
else {
if (rank==0) printf("Test failed!: mass conservation for D3Q7 \n");
}
}
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();

2
tests/lbpm_color_simulator.cpp
View File

@@ -58,7 +58,7 @@ int main(int argc, char **argv)
ColorModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
ColorModel.Run();
//ColorModel.WriteDebug();
ColorModel.WriteDebug();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_color_simulator",1);