diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 3047bd0a..9c2f4529 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -245,8 +245,8 @@ void SubPhase::Basic(){ force_mag = 1.0; } double saturation=gwb.V/(gwb.V + gnb.V); - double water_flow_rate=gwb.V*sqrt(gwb.Px*gwb.Px + gwb.Py*gwb.Py + gwb.Pz*gwb.Pz)/gwb.M / Dm->Volume; - double not_water_flow_rate=gnb.V*sqrt(gnb.Px*gnb.Px + gnb.Py*gnb.Py + gnb.Pz*gnb.Pz)/gnb.M/ Dm->Volume; + double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume; + double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; double total_flow_rate = water_flow_rate + not_water_flow_rate; double fractional_flow= water_flow_rate / total_flow_rate; diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu index ae96efb3..2df4db6c 100644 --- a/gpu/D3Q19.cu +++ b/gpu/D3Q19.cu @@ -2499,9 +2499,8 @@ extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double // Allocate memory to store the sums double din; - double *sum; + double sum[1]; double *dvcsum; - cudaMallocHost((void **)&sum,sizeof(double)); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2544,9 +2543,8 @@ extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, dou // Allocate memory to store the sums double din; - double *sum; + double sum[1]; double *dvcsum; - cudaMallocHost((void **)&sum,sizeof(double)); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2595,7 +2593,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, doub dvc_D3Q19_Flux_BC_Z<<>>(disteven, distodd, flux, dvcsum, Nx, Ny, Nz, outlet); // Now read the total flux - cudaMemcpy(&sum[0],&dvcsum[0],sizeof(double),cudaMemcpyDeviceToHost); + cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost); // free the memory needed for reduction diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index d6511b7b..7861753b 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -654,8 +654,8 @@ void ScaLBL_ColorModel::Run(){ force_mag = 1.0; } double current_saturation = volB/(volA+volB); - double flow_rate_A = volA*sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); - double flow_rate_B = volB*sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); + double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z); + double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z); double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha); if ( morph_timesteps > morph_interval ){ diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index d4fb775c..f88f2de4 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -23,7 +23,7 @@ void ScaLBL_MRTModel::ReadParams(string filename){ tau = 1.0; timestepMax = 100000; - tolerance = 0.01; + tolerance = 1.0e-8; Fx = Fy = 0.0; Fz = 1.0e-5; @@ -203,12 +203,12 @@ void ScaLBL_MRTModel::Run(){ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); MPI_Barrier(comm); starttime = MPI_Wtime(); - if (rank==0) printf("Beginning AA timesteps...\n"); + if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax); if (rank==0) printf("********************************************************\n"); timestep=0; double error = 1.0; double flow_rate_previous = 0.0; - while (timestep < timestepMax && error < tolerance) { + while (timestep < timestepMax && error > tolerance) { //************************************************************************/ timestep++; ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9ae38153..9737c7ad 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -41,6 +41,7 @@ CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_ ADD_LBPM_TEST( pmmc_cylinder ) ADD_LBPM_TEST( TestTorus ) ADD_LBPM_TEST( TestTorusEvolve ) +ADD_LBPM_TEST( TestTopo3D ) ADD_LBPM_TEST( TestFluxBC ) ADD_LBPM_TEST( TestMap ) #ADD_LBPM_TEST( TestMRT ) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp new file mode 100644 index 00000000..4b8cd17f --- /dev/null +++ b/tests/TestTopo3D.cpp @@ -0,0 +1,234 @@ +// Sequential blob analysis +// Reads parallel simulation data and performs connectivity analysis +// and averaging on a blob-by-blob basis +// James E. McClure 2014 + +#include +#include +#include "common/Communication.h" +#include "analysis/analysis.h" +#include "analysis/Minkowski.h" +#include "IO/MeshDatabase.h" + +std::shared_ptr loadInputs( int nprocs ) +{ + //auto db = std::make_shared( "Domain.in" ); + auto db = std::make_shared(); + db->putScalar( "BC", 0 ); + db->putVector( "nproc", { 1, 1, 1 } ); + db->putVector( "n", { 100, 100, 100 } ); + db->putScalar( "nspheres", 1 ); + db->putVector( "L", { 1, 1, 1 } ); + return db; +} + + +int main(int argc, char **argv) +{ + // Initialize MPI + int rank, nprocs; + MPI_Init(&argc,&argv); + MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + if ( rank==0 ) { + printf("-----------------------------------------------------------\n"); + printf("Unit test 3D topologies \n"); + printf("-----------------------------------------------------------\n"); + } + + //....................................................................... + // Reading the domain information file + //....................................................................... + int i,j,k,n; + + // Load inputs + auto db = loadInputs( nprocs ); + int Nx = db->getVector( "n" )[0]; + int Ny = db->getVector( "n" )[1]; + int Nz = db->getVector( "n" )[2]; + int nprocx = db->getVector( "nproc" )[0]; + int nprocy = db->getVector( "nproc" )[1]; + int nprocz = db->getVector( "nproc" )[2]; + + if (rank==0){ + printf("********************************************************\n"); + printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz); + printf("********************************************************\n"); + } + + // Get the rank info + std::shared_ptr Dm(new Domain(db,comm)); + + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + //....................................................................... + for ( k=1;kid[n] = 1; + } + } + } + //....................................................................... + Dm->CommInit(); // Initialize communications for domains + //....................................................................... + + // Create visualization structure + std::vector visData; + fillHalo fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);; + + IO::initialize("","silo","false"); + // Create the MeshDataStruct + visData.resize(1); + visData[0].meshName = "domain"; + visData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); + auto PhaseVar = std::make_shared(); + PhaseVar->name = "phase"; + PhaseVar->type = IO::VariableType::VolumeVariable; + PhaseVar->dim = 1; + PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(PhaseVar); + + //....................................................................... + // Assign the phase ID field based and the signed distance + //....................................................................... + double R1,R2,R; + double CX,CY,CZ; //CY1,CY2; + CX=Nx*nprocx*0.5; + CY=Ny*nprocy*0.5; + CZ=Nz*nprocz*0.5; + R1 = (Nx-2)*nprocx*0.3; // middle radius + R2 = (Nx-2)*nprocx*0.1; // donut thickness + R = 0.4*nprocx*(Nx-2); + + Minkowski Object(Dm); + + int timestep = 0; + double x,y,z; + + // partial torus + timestep += 1; + for ( k=1;kiproc()*(Nx-2)+i - CX - 0.1; + y = Dm->jproc()*(Ny-2)+j - CY - 0.1; + z = Dm->kproc()*(Nz-2)+k - CZ -0.1; + + //.............................................................................. + if (x <= 0 || y<=0) { + // Single torus + Object.distance(i,j,k) = R2 - sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z); + } + else { + double d1 = R2-sqrt(x*x +(y-R1)*(y-R1) + z*z); + double d2 = R2-sqrt((x-R1)*(x-R1)+y*y + z*z); + Object.distance(i,j,k) = max(d1,d2); + } + + if (Object.distance(i,j,k) > 0.0){ + Dm->id[n] = 2; + Object.id(i,j,k) = 2; + } + else{ + Dm->id[n] = 1; + Object.id(i,j,k) = 1; + } + } + } + } + + ASSERT(visData[0].vars[0]->name=="phase"); + Array& PhaseData = visData[0].vars[0]->data; + fillData.copy(Object.distance,PhaseData); + IO::writeData( timestep, visData, comm ); + + //spherical shell + timestep += 1; + for ( k=1;kiproc()*(Nx-2)+i - CX - 0.1; + y = Dm->jproc()*(Ny-2)+j - CY - 0.1; + z = Dm->kproc()*(Nz-2)+k - CZ - 0.1; + + //.............................................................................. + // Single torus + double d1 = sqrt(x*x+y*y+z*z)-(R1-R2); + double d2 = R-sqrt(x*x+y*y+z*z); + Object.distance(i,j,k) = min(d1,d2); + + if (Object.distance(i,j,k) > 0.0){ + Dm->id[n] = 2; + Object.id(i,j,k) = 2; + } + else{ + Dm->id[n] = 1; + Object.id(i,j,k) = 1; + } + } + } + } + ASSERT(visData[0].vars[0]->name=="phase"); + PhaseData = visData[0].vars[0]->data; + fillData.copy(Object.distance,PhaseData); + IO::writeData( timestep, visData, comm ); + + // bowl + timestep += 1; + for ( k=1;kiproc()*(Nx-2)+i - CX - 0.1; + y = Dm->jproc()*(Ny-2)+j - CY - 0.1; + z = Dm->kproc()*(Nz-2)+k - CZ - 0.1; + + //.............................................................................. + // Bowl + if (z > 0 ){ + Object.distance(i,j,k) = R2-sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z); + } + else + { + double d1 = sqrt(x*x+y*y+z*z)-(R1-R2); + double d2 = R-sqrt(x*x+y*y+z*z); + Object.distance(i,j,k) = min(d1,d2); + } + if (Object.distance(i,j,k) > 0.0){ + Dm->id[n] = 2; + Object.id(i,j,k) = 2; + } + else{ + Dm->id[n] = 1; + Object.id(i,j,k) = 1; + } + } + } + } + + ASSERT(visData[0].vars[0]->name=="phase"); + PhaseData = visData[0].vars[0]->data; + fillData.copy(Object.distance,PhaseData); + IO::writeData( timestep, visData, comm ); + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + MPI_Barrier(comm); + MPI_Finalize(); + return 0; +} +