partial merge
This commit is contained in:
James McClure
@@ -8,7 +8,7 @@
|
||||
#include <fstream>
|
||||
|
||||
#include "common/ScaLBL.h"
|
||||
#include "common/MPI.h"
|
||||
#include "common/MPI_Helpers.h"
|
||||
#include "models/ColorModel.h"
|
||||
|
||||
inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
|
||||
@@ -41,7 +41,7 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
|
||||
int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
|
||||
int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
|
||||
// Initialize phase position field for parallel bubble test
|
||||
if (jglobal < 40){
|
||||
if (kglobal < 40){
|
||||
ColorModel.Mask->id[n] = 0;
|
||||
}
|
||||
else if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
|
||||
@@ -72,6 +72,7 @@ int main(int argc, char **argv)
|
||||
int rank = comm.getRank();
|
||||
int nprocs = comm.getSize();
|
||||
// parallel domain size (# of sub-domains)
|
||||
int CleanCheck = 0;
|
||||
|
||||
if (rank == 0){
|
||||
printf("********************************************************\n");
|
||||
@@ -83,68 +84,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];
|
||||
@@ -152,60 +153,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 (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-8 || fabs(total_mass_B_1-total_mass_B_0) > 1.0e-8) 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
|
||||
@@ -255,14 +257,14 @@ 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");
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
Utilities::shutdown();
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user