From 081342d8bdf8110170fabdb866ee0586279c5598 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 30 Nov 2018 16:45:59 -0500 Subject: [PATCH] working on mass conservation unit test --- tests/TestMassConservationD3Q7.cpp | 220 +++++++++++++++-------------- 1 file changed, 115 insertions(+), 105 deletions(-) diff --git a/tests/TestMassConservationD3Q7.cpp b/tests/TestMassConservationD3Q7.cpp index a08342a1..e37b93fc 100644 --- a/tests/TestMassConservationD3Q7.cpp +++ b/tests/TestMassConservationD3Q7.cpp @@ -8,9 +8,61 @@ #include #include "common/ScaLBL.h" -#include "common/Communication.h" -#include "analysis/TwoPhase.h" #include "common/MPI_Helpers.h" +#include "models/ColorModel.h" + +inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){ + // initialize a bubble + int i,j,k,n; + int rank = ColorModel.Mask->rank(); + int nprocx = ColorModel.Mask->nprocx(); + int nprocy = ColorModel.Mask->nprocy(); + int nprocz = ColorModel.Mask->nprocz(); + int Nx = ColorModel.Mask->Nx; + int Ny = ColorModel.Mask->Ny; + int Nz = ColorModel.Mask->Nz; + if (rank == 0) cout << "Setting up bubble..." << endl; + for (k=0;kSDs(i,j,k) = 100.f; + } + } + } + // Initialize the bubble + int count_in_bubble=0; + int count_out_bubble=0; + for (k=0;kiproc(); + 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){ + ColorModel.Mask->id[n] = 0; + } + else if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx) + +(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy) + +(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){ + ColorModel.Mask->id[n] = 2; + ColorModel.Mask->id[n] = 2; + count_in_bubble++; + } + else{ + ColorModel.Mask->id[n]=1; + ColorModel.Mask->id[n]=1; + count_out_bubble++; + } + ColorModel.id[n] = ColorModel.Mask->id[n]; + } + } + } + printf("Count in %i, out %i\n",count_in_bubble,count_out_bubble); + // initialize the phase indicator field +} int main(int argc, char **argv) { @@ -42,130 +94,88 @@ int main(int argc, char **argv) printf("********************************************************\n"); printf("Running Unit Test for D3Q7 Mass Conservation \n"); printf("********************************************************\n"); - } - - double das=0.0; - double dbs=1.0; - double beta=0.95; - bool pBC=false; - int i,j,k,n; - int Nx,Ny,Nz,N; - Nx=Nz=Ny=100; - N = Nx*Ny*Nz; - int dist_mem_size = N*sizeof(double); - - double *DenOriginal, *DenFinal; - DenOriginal = new double [2*N]; - DenFinal = new double [2*N]; - - double *Vel; - Vel = new double [3*N]; - char *id; - id = new char [N]; - - for (k=0;k80) id[n] = 0; - else if (i<50){ - id[n]=1; - Vel[n]=0.1; - Vel[N+n]=0.1; - Vel[2*N+n]=0.1; - } - else { - id[n]=2; - Vel[n]=0.1; - Vel[N+n]=0.1; - Vel[2*N+n]=0.1; - } - } + if ( argc < 2 ) { + std::cerr << "Invalid number of arguments, no input file specified\n"; + return -1; } } + { + 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; + Np = CM.Np; + N = Nx*Ny*Nz; + int dist_mem_size = N*sizeof(double); + double *DenOriginal, *DenFinal; + DenOriginal = new double [2*Np]; + DenFinal = new double [2*Np]; + + //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(); - - //......................device distributions................................. - double *A_even,*A_odd,*B_even,*B_odd; - //........................................................................... - ScaLBL_AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size); // Allocate device memory - ScaLBL_AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size); // Allocate device memory - ScaLBL_AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size); // Allocate device memory - ScaLBL_AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size); // Allocate device memory - //........................................................................... - double *Phi,*Den; - double *ColorGrad, *Velocity; - //........................................................................... - ScaLBL_AllocateDeviceMemory((void **) &Phi, dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size); - //...........device phase ID................................................. - if (rank==0) printf ("Copying phase ID to device \n"); - char *ID; - ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory - // Copy to the device - ScaLBL_CopyToDevice(ID, id, N); - ScaLBL_CopyToDevice(Velocity, Vel, 3*N*sizeof(double)); - //........................................................................... - - ScaLBL_Color_Init(ID, Den, Phi, das, dbs, Nx, Ny, Nz); - ScaLBL_CopyToHost(DenOriginal,Den,2*N*sizeof(double)); + // Run the odd timestep + ScaLBL_CopyToHost(DenOriginal,CM.Den,2*Np*sizeof(double)); - //...................................................................... - ScaLBL_D3Q7_Init(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz); - ScaLBL_D3Q7_Init(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz); - ScaLBL_ComputePhaseField(ID, Phi, Den, N); - ScaLBL_D3Q19_ColorGradient(ID,Phi,ColorGrad,Nx,Ny,Nz); - //.................................................................................. + 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); - //************************************************************************* - // Carry out the density streaming step for mass transport - //************************************************************************* - ScaLBL_D3Q7_ColorCollideMass(ID, A_even, A_odd, B_even, B_odd, Den, Phi, - ColorGrad, Velocity, beta, N, pBC); - //************************************************************************* - - //.................................................................................. - ScaLBL_D3Q7_Density(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz); - ScaLBL_D3Q7_Density(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz); - //.................................................................................. - ScaLBL_ComputePhaseField(ID, Phi, Den, N); - ScaLBL_D3Q19_ColorGradient(ID,Phi,ColorGrad,Nx,Ny,Nz); - //.................................................................................. + FILE *OUTFILE; + OUTFILE = fopen("phi-out.raw","wb"); + fwrite(CM.Phi,8,N,OUTFILE); + fclose(OUTFILE); // Compare and make sure mass is conserved at every lattice site + double *Error; + Error = new double [N]; bool CleanCheck = true; double original,final; - ScaLBL_CopyToHost(DenFinal,Den,2*N*sizeof(double)); + ScaLBL_CopyToHost(DenFinal,CM.Den,2*Np*sizeof(double)); for (k=0;k 1e-15){ - final = DenFinal[n]; - original = DenOriginal[n]; - if (id[n] == 0) printf("Solid phase! \n"); - if (id[n] == 1) printf("Wetting phase! \n"); - if (id[n] == 2) printf("Non-wetting phase! \n"); + + int idx = CM.Map(i,j,k); + if (idx < Np && !(idx<0)){ + printf("idx=%i\n",idx); + if (fabs(DenFinal[idx] - DenOriginal[idx]) > 1e-15){ + final = DenFinal[idx]; + original = DenOriginal[idx]; + //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; } - if (fabs(DenFinal[N+n] - DenOriginal[N+n]) > 1e-15){ - if (id[n] == 0) printf("Solid phase! \n"); - if (id[n] == 1) printf("Wetting phase! \n"); - if (id[n] == 2) printf("Non-wetting phase! \n"); - final = DenFinal[N+n]; - original = DenOriginal[N+n]; + 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"); + final = DenFinal[Np+idx]; + original = DenOriginal[Np+idx]; printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final); CleanCheck=false; } + } } } } - if (rank==0) printf("Checking that the correct velocity is retained \n"); +/* if (rank==0) printf("Checking that the correct velocity is retained \n"); // Swap convention is observed -- velocity is negative double *Aeven,*Aodd,*Beven,*Bodd; Aeven = new double[4*N]; @@ -213,7 +223,7 @@ int main(int argc, char **argv) } } } - +*/ if (CleanCheck){ if (rank==0) printf("Test passed: mass conservation for D3Q7 \n"); } @@ -221,7 +231,7 @@ int main(int argc, char **argv) if (rank==0) printf("Test failed!: mass conservation for D3Q7 \n"); } - +} // **************************************************** MPI_Barrier(comm); MPI_Finalize();