From 84f24d8fc5e4e3c5fae76248f7688968ba26c862 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 20 May 2019 09:46:55 -0600 Subject: [PATCH 01/11] debugging flux bc --- gpu/D3Q19.cu | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu index ae96efb3..afbad77c 100644 --- a/gpu/D3Q19.cu +++ b/gpu/D3Q19.cu @@ -2501,7 +2501,8 @@ extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double double din; double *sum; double *dvcsum; - cudaMallocHost((void **)&sum,sizeof(double)); + sum = new double [count]; + //cudaMallocHost((void **)&sum,sizeof(double)*count); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2519,7 +2520,7 @@ extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double } // Now read the total flux - cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost); + cudaMemcpy(&sum,dvcsum,sizeof(double)*count,cudaMemcpyDeviceToHost); din=sum[0]; err = cudaGetLastError(); if (cudaSuccess != err){ @@ -2546,7 +2547,8 @@ extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, dou double din; double *sum; double *dvcsum; - cudaMallocHost((void **)&sum,sizeof(double)); + sum = new double [count]; + //cudaMallocHost((void **)&sum,sizeof(double)*count); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2562,7 +2564,7 @@ extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, dou printf("CUDA error in ScaLBL_D3Q19_AAodd_Flux_BC_z (kernel): %s \n",cudaGetErrorString(err)); } // Now read the total flux - cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost); + cudaMemcpy(&sum,dvcsum,sizeof(double)*count,cudaMemcpyDeviceToHost); din=sum[0]; err = cudaGetLastError(); if (cudaSuccess != err){ From 79c2d2704c5da8363a3dbc28a60e8319bdbaaa96 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 20 May 2019 09:57:54 -0600 Subject: [PATCH 02/11] relperm aligned with flow direction --- analysis/SubPhase.cpp | 4 ++-- models/ColorModel.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) 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/models/ColorModel.cpp b/models/ColorModel.cpp index 3d64a0fd..fdff6dba 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -655,8 +655,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 ){ From 73c80fcbc78d28f56c37e080b43137689d632040 Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 20 May 2019 20:37:50 -0600 Subject: [PATCH 03/11] reverting flux bc changes --- gpu/D3Q19.cu | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu index afbad77c..2df4db6c 100644 --- a/gpu/D3Q19.cu +++ b/gpu/D3Q19.cu @@ -2499,10 +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; - sum = new double [count]; - //cudaMallocHost((void **)&sum,sizeof(double)*count); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2520,7 +2518,7 @@ extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double } // Now read the total flux - cudaMemcpy(&sum,dvcsum,sizeof(double)*count,cudaMemcpyDeviceToHost); + cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost); din=sum[0]; err = cudaGetLastError(); if (cudaSuccess != err){ @@ -2545,10 +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; - sum = new double [count]; - //cudaMallocHost((void **)&sum,sizeof(double)*count); cudaMalloc((void **)&dvcsum,sizeof(double)*count); cudaMemset(dvcsum,0,sizeof(double)*count); int sharedBytes = 512*sizeof(double); @@ -2564,7 +2560,7 @@ extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, dou printf("CUDA error in ScaLBL_D3Q19_AAodd_Flux_BC_z (kernel): %s \n",cudaGetErrorString(err)); } // Now read the total flux - cudaMemcpy(&sum,dvcsum,sizeof(double)*count,cudaMemcpyDeviceToHost); + cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost); din=sum[0]; err = cudaGetLastError(); if (cudaSuccess != err){ @@ -2597,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 From 176acd3421103d9240064d793ac1f62d2bf952ef Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 16:50:56 -0600 Subject: [PATCH 04/11] add geometric test for 3D topo --- tests/CMakeLists.txt | 1 + tests/TestTopo3D.cpp | 232 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 233 insertions(+) create mode 100644 tests/TestTopo3D.cpp 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..c0fb1e3a --- /dev/null +++ b/tests/TestTopo3D.cpp @@ -0,0 +1,232 @@ +// 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){ + double d1 = R2-sqrt(x*x +(y-R1)*(y-R1)); + double d2 = R2-sqrt((x-R1)*(x-R1)+y*y); + Object.distance(i,j,k) = max(d1,d2); + } + else{ + // Single torus + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } + + 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; + 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"); + Array& 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 ){ + double d1 = sqrt(x*x+y*y+z*z)-R1; + double d2 = R-sqrt(x*x+y*y+z*z); + Object.distance(i,j,k) = min(d1,d2); + } + else{ + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } + + 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 ); + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + MPI_Barrier(comm); + MPI_Finalize(); + return 0; +} + From a3afb35114e9f423128701fda03b6d721e07e5c7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 16:55:27 -0600 Subject: [PATCH 05/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 388 ++++++++++++++++++++++--------------------- 1 file changed, 195 insertions(+), 193 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index c0fb1e3a..3748a408 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -12,221 +12,223 @@ 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; + //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 + // 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; - } + if ( rank==0 ) { + printf("-----------------------------------------------------------\n"); + printf("Unit test 3D topologies \n"); + printf("-----------------------------------------------------------\n"); } - } - //....................................................................... - 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); + //....................................................................... + // Reading the domain information file + //....................................................................... + int i,j,k,n; - int timestep = 0; - double x,y,z; + // 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]; - // 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; + // Get the rank info + std::shared_ptr Dm(new Domain(db,comm)); - //.............................................................................. - if (x > 0 && y>0){ - double d1 = R2-sqrt(x*x +(y-R1)*(y-R1)); - double d2 = R2-sqrt((x-R1)*(x-R1)+y*y); - Object.distance(i,j,k) = max(d1,d2); - } - else{ - // Single torus - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; - } + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + //....................................................................... + for ( k=1;kid[n] = 1; + } + } + } + //....................................................................... + Dm->CommInit(); // Initialize communications for domains + //....................................................................... - 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; - } - } - } + // 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);; - 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;k( 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); - // global position relative to center - x = Dm->iproc()*(Nx-2)+i - CX - 0.1; - y = Dm->jproc()*(Ny-2)+j - CY - 0.1; - z = Dm->kproc()*(Nz-2)+k - CZ - 0.1; + //....................................................................... + // 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); - //.............................................................................. - // Single torus - double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = R-sqrt(x*x+y*y+z*z); - Object.distance(i,j,k) = min(d1,d2); + Minkowski Object(Dm); - 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; - } - } - } + int timestep = 0; + double x,y,z; - ASSERT(visData[0].vars[0]->name=="phase"); - Array& 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; + // global position relative to center + x = Dm->iproc()*(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 ){ - double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = R-sqrt(x*x+y*y+z*z); - Object.distance(i,j,k) = min(d1,d2); - } - else{ - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; - } + //.............................................................................. + if (x > 0 && y>0){ + double d1 = R2-sqrt(x*x +(y-R1)*(y-R1)); + double d2 = R2-sqrt((x-R1)*(x-R1)+y*y); + Object.distance(i,j,k) = max(d1,d2); + } + else{ + // Single torus + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } - 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; - } - } - } + 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 ); + ASSERT(visData[0].vars[0]->name=="phase"); + Array& 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; + //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; + 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"); + Array& 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 ){ + double d1 = sqrt(x*x+y*y+z*z)-R1; + double d2 = R-sqrt(x*x+y*y+z*z); + Object.distance(i,j,k) = min(d1,d2); + } + else{ + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } + + 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 ); + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + MPI_Barrier(comm); + MPI_Finalize(); + return 0; } From 5042b4f0cd41ca5c9f1fdb0d3c37653370961cca Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 16:56:51 -0600 Subject: [PATCH 06/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index 3748a408..4c724637 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -182,7 +182,7 @@ int main(int argc, char **argv) } } ASSERT(visData[0].vars[0]->name=="phase"); - Array& PhaseData = visData[0].vars[0]->data; + PhaseData = visData[0].vars[0]->data; fillData.copy(Object.distance,PhaseData); IO::writeData( timestep, visData, comm ); @@ -222,7 +222,7 @@ int main(int argc, char **argv) } ASSERT(visData[0].vars[0]->name=="phase"); - Array& PhaseData = visData[0].vars[0]->data; + PhaseData = visData[0].vars[0]->data; fillData.copy(Object.distance,PhaseData); IO::writeData( timestep, visData, comm ); From 9161b2caa65d1151533ed52475228623aba4ac35 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 17:09:01 -0600 Subject: [PATCH 07/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index 4c724637..8bdcd71f 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -126,8 +126,8 @@ int main(int argc, char **argv) //.............................................................................. if (x > 0 && y>0){ - double d1 = R2-sqrt(x*x +(y-R1)*(y-R1)); - double d2 = R2-sqrt((x-R1)*(x-R1)+y*y); + 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); } else{ @@ -167,7 +167,7 @@ int main(int argc, char **argv) //.............................................................................. // Single torus double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = R-sqrt(x*x+y*y+z*z); + double d2 = (R1-R2)-sqrt(x*x+y*y+z*z); Object.distance(i,j,k) = min(d1,d2); if (Object.distance(i,j,k) > 0.0){ @@ -200,9 +200,9 @@ int main(int argc, char **argv) //.............................................................................. // Bowl - if (z <0 ){ + if (z < 0 ){ double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = R-sqrt(x*x+y*y+z*z); + double d2 = (R1-R2)-sqrt(x*x+y*y+z*z); Object.distance(i,j,k) = min(d1,d2); } else{ From d1cec4fb9b5c67f3532c721249b49931dfe1183e Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 17:14:16 -0600 Subject: [PATCH 08/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index 8bdcd71f..0e81815e 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -125,7 +125,7 @@ int main(int argc, char **argv) z = Dm->kproc()*(Nz-2)+k - CZ -0.1; //.............................................................................. - if (x > 0 && y>0){ + if (x >= 0 && y>=0){ 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); @@ -166,8 +166,8 @@ int main(int argc, char **argv) //.............................................................................. // Single torus - double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = (R1-R2)-sqrt(x*x+y*y+z*z); + 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){ @@ -200,9 +200,9 @@ int main(int argc, char **argv) //.............................................................................. // Bowl - if (z < 0 ){ - double d1 = sqrt(x*x+y*y+z*z)-R1; - double d2 = (R1-R2)-sqrt(x*x+y*y+z*z); + if (z <= 0 ){ + 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); } else{ From 06ecb0db4e492624e781859c969a725044dee39f Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 17:21:25 -0600 Subject: [PATCH 09/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index 0e81815e..e7676a67 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -125,15 +125,15 @@ int main(int argc, char **argv) z = Dm->kproc()*(Nz-2)+k - CZ -0.1; //.............................................................................. - if (x >= 0 && y>=0){ + if (x <= 0 || y<=0) { + // Single torus + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } + 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); } - else{ - // Single torus - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; - } if (Object.distance(i,j,k) > 0.0){ Dm->id[n] = 2; @@ -200,15 +200,15 @@ int main(int argc, char **argv) //.............................................................................. // Bowl - if (z <= 0 ){ + if (z > 0 ){ + Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + } + 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); } - else{ - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; - } - if (Object.distance(i,j,k) > 0.0){ Dm->id[n] = 2; Object.id(i,j,k) = 2; From f360eb1ecb9c3184ea63757fc5c65003abaecf20 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 22 May 2019 17:33:10 -0600 Subject: [PATCH 10/11] add geometric test for 3D topo --- tests/TestTopo3D.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/TestTopo3D.cpp b/tests/TestTopo3D.cpp index e7676a67..4b8cd17f 100644 --- a/tests/TestTopo3D.cpp +++ b/tests/TestTopo3D.cpp @@ -127,7 +127,7 @@ int main(int argc, char **argv) //.............................................................................. if (x <= 0 || y<=0) { // Single torus - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + 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); @@ -201,7 +201,7 @@ int main(int argc, char **argv) //.............................................................................. // Bowl if (z > 0 ){ - Object.distance(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + Object.distance(i,j,k) = R2-sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z); } else { From 8ad659ff17f31efa810162cc200e65adb8f197d9 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 5 Jun 2019 14:57:06 -0400 Subject: [PATCH 11/11] fix bug in permeability tolerance --- models/MRTModel.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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