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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2018-12-02 07:48:44 -05:00
commit 807e476bea
4 changed files with 213 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
models/ColorModel.cpp
View File

@ -184,15 +184,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n");
}
if (rank==0){
printf("Components labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
printf(" label=%i, affinity=%f\n",int(VALUE),AFFINITY);
}
}
int label_count[NLABELS];
int label_count_global[NLABELS];
// Assign the labels
for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
@ -200,9 +197,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
VALUE=id[n];
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("rank=%i, idx=%i, value=%i, %i, \n",rank(),idx, VALUE,LabelList[idx]);
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
label_count[idx]++;
idx = NLABELS;
Mask->id[n] = 0; // set mask to zero since this is an immobile component
}
@ -216,6 +214,19 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
}
// Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
MPI_Allreduce(&label_count[0],&label_count_global[0],NLABELS,MPI_INT,MPI_SUM,Dm->Comm);
if (rank==0){
printf("Components labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
double volume_fraction = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
printf(" label=%i, affinity=%f, volume fraction==%f\n",int(VALUE),AFFINITY,volume_fraction);
}
}
}
@ -592,8 +603,8 @@ void ScaLBL_ColorModel::Run(){
if (rank==0){
printf("** WRITE STEADY POINT *** ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
volA /= double(Nx*Ny*Nz*nprocs);
volB /= double(Nx*Ny*Nz*nprocs);
volA /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
volB /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
FILE * kr_log_file = fopen("relperm.csv","a");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g ",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz);
fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z);
@ -660,10 +671,17 @@ void ScaLBL_ColorModel::Run(){
if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum
if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta);
}
//MORPH_ADAPT = false;
if (volB/(volA + volB) > TARGET_SATURATION){
MORPH_ADAPT = false;
TARGET_SATURATION = target_saturation[target_saturation_index++];
if (morph_delta < 0.f){
if (volB/(volA + volB) > TARGET_SATURATION){
MORPH_ADAPT = false;
TARGET_SATURATION = target_saturation[target_saturation_index++];
}
}
else{
if (volB/(volA + volB) < TARGET_SATURATION){
MORPH_ADAPT = false;
TARGET_SATURATION = target_saturation[target_saturation_index++];
}
}
MPI_Barrier(comm);
morph_timesteps = 0;

2
tests/CMakeLists.txt
View File

@ -56,13 +56,13 @@ ADD_LBPM_TEST( TestMomentsD3Q19 )
ADD_LBPM_TEST( TestInterfaceSpeed ../example/Bubble/input.db)
ADD_LBPM_TEST( test_dcel_minkowski )
ADD_LBPM_TEST( test_dcel_tri_normal )
ADD_LBPM_TEST( TestMassConservationD3Q7 ../example/Bubble/input.db)
#ADD_LBPM_TEST_1_2_4( TestTwoPhase )
ADD_LBPM_TEST_1_2_4( TestBlobIdentify )
#ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )
#ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 )
ADD_LBPM_TEST_PARALLEL( TestSegDist 8 )
ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 )
#ADD_LBPM_TEST_PARALLEL( TestMassConservationD3Q7 1 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST_1_2_4( testUtilities )
ADD_LBPM_TEST( TestWriter )

288
tests/TestMassConservationD3Q7.cpp
View File

@ -8,9 +8,61 @@
#include <fstream>
#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;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
ColorModel.Averages->SDs(i,j,k) = 100.f;
}
}
}
// Initialize the bubble
int count_in_bubble=0;
int count_out_bubble=0;
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;
int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
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,129 +94,139 @@ int main(int argc, char **argv)
printf("********************************************************\n");
printf("Running Unit Test for D3Q7 Mass Conservation \n");
printf("********************************************************\n");
if ( argc < 2 ) {
std::cerr << "Invalid number of arguments, no input file specified\n";
return -1;
}
}
double das=0.0;
double dbs=1.0;
double beta=0.95;
bool pBC=false;
{
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;
Nx=Nz=Ny=100;
int Nx,Ny,Nz,N,Np;
Nx = CM.Nx;
Ny = CM.Ny;
Nz = CM.Nz;
N = Nx*Ny*Nz;
int dist_mem_size = N*sizeof(double);
//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();
Np = CM.Np;
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;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
if (k<10) id[n] = 0;
else if (k>80) 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;
}
}
}
}
DenOriginal = new double [2*Np];
DenFinal = new double [2*Np];
//......................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));
//......................................................................
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);
//..................................................................................
//*************************************************************************
// 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);
//..................................................................................
// 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
double *Error;
Error = new double [N];
double *A_q, *B_q;
A_q = new double [7*Np];
B_q = new double [7*Np];
bool CleanCheck = true;
double original,final;
ScaLBL_CopyToHost(DenFinal,Den,2*N*sizeof(double));
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
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;
if (fabs(DenFinal[n] - DenOriginal[n]) > 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");
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];
printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
CleanCheck=false;
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];
total_mass_A_0 += original;
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];
}
}
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;
}*/
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;
}*/
}
}
}
}
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;
/*
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
double *Aeven,*Aodd,*Beven,*Bodd;
@ -213,7 +275,7 @@ int main(int argc, char **argv)
}
}
}
*/
if (CleanCheck){
if (rank==0) printf("Test passed: mass conservation for D3Q7 \n");
}
@ -221,7 +283,7 @@ int main(int argc, char **argv)
if (rank==0) printf("Test failed!: mass conservation for D3Q7 \n");
}
}
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();

4
tests/lbpm_serial_decomp.cpp
View File

@ -168,6 +168,10 @@ int main(int argc, char **argv)
int64_t y = yStart + jp*ny + j-1;
// int64_t z = zStart + kp*nz + k-1;
int64_t z = zStart + kp*nz + k-1 - z_transition_size;
if (x<xStart) x=xStart;
if (!(x<Nx)) x=Nx-1;
if (y<yStart) y=yStart;
if (!(y<Ny)) y=Ny-1;
if (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;