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; }