diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 9c2f4529..a9d6f910 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -231,7 +231,20 @@ void SubPhase::Basic(){ count_n=sumReduce( Dm->Comm, count_n); gwb.p=sumReduce( Dm->Comm, wb.p) / count_w; gnb.p=sumReduce( Dm->Comm, nb.p) / count_n; - + + // check for NaN + bool err=false; + if (gwb.V != gwb.V) err=true; + if (gnb.V != gnb.V) err=true; + if (gwb.p != gwb.p) err=true; + if (gnb.p != gnb.p) err=true; + if (gwb.Px != gwb.Px) err=true; + if (gwb.Py != gwb.Py) err=true; + if (gwb.Pz != gwb.Pz) err=true; + if (gnb.Px != gnb.Px) err=true; + if (gnb.Py != gnb.Py) err=true; + if (gnb.Pz != gnb.Pz) err=true; + if (Dm->rank() == 0){ double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); double dir_x = Fx/force_mag; @@ -257,7 +270,11 @@ void SubPhase::Basic(){ fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p); fflush(TIMELOG); } - + if (err==true){ + // exception if simulation produceds NaN + printf("SubPhase.cpp: NaN encountered, may need to check simulation parameters \n"); + } + ASSERT(err==false); } inline void InterfaceTransportMeasures( double beta, double rA, double rB, double nA, double nB, @@ -492,6 +509,10 @@ void SubPhase::Full(){ double vol_wc_bulk = 0.0; double vol_nd_bulk = 0.0; double vol_wd_bulk = 0.0; + double count_wc = 0.0; + double count_nc = 0.0; + double count_wd = 0.0; + double count_nd = 0.0; for (k=kmin; k 0.0){ + if (morph_n->label(i,j,k) > 0 ){ + count_nd += 1.0; + nd.p += Pressure(n); + } + else{ + count_nc += 1.0; + nc.p += Pressure(n); + } + } + else{ + // water region + if (morph_w->label(i,j,k) > 0 ){ + count_wd += 1.0; + wd.p += Pressure(n); + } + else{ + count_wc += 1.0; + wc.p += Pressure(n); + } + } + if ( phi > 0.0){ if (morph_n->label(i,j,k) > 0 ){ vol_nd_bulk += 1.0; nd.M += nA*rho_n; @@ -522,7 +564,6 @@ void SubPhase::Full(){ nd.Py += nA*rho_n*uy; nd.Pz += nA*rho_n*uz; nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); - nd.p += Pressure(n); } else{ vol_nc_bulk += 1.0; @@ -531,7 +572,6 @@ void SubPhase::Full(){ nc.Py += nA*rho_n*uy; nc.Pz += nA*rho_n*uz; nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); - nc.p += Pressure(n); } } else{ @@ -543,7 +583,6 @@ void SubPhase::Full(){ wd.Py += nB*rho_w*uy; wd.Pz += nB*rho_w*uz; wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); - wd.p += Pressure(n); } else{ vol_wc_bulk += 1.0; @@ -552,40 +591,44 @@ void SubPhase::Full(){ wc.Py += nB*rho_w*uy; wc.Pz += nB*rho_w*uz; wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); - wc.p += Pressure(n); } } } } } } + count_wc=sumReduce( Dm->Comm, count_wc); + count_nc=sumReduce( Dm->Comm, count_nc); + count_wd=sumReduce( Dm->Comm, count_wd); + count_nd=sumReduce( Dm->Comm, count_nd); + gnd.p=sumReduce( Dm->Comm, nd.p) / count_nd; + gwd.p=sumReduce( Dm->Comm, wd.p) / count_wd; + gnc.p=sumReduce( Dm->Comm, nc.p) / count_nc; + gwc.p=sumReduce( Dm->Comm, wc.p) / count_wc; + gnd.M=sumReduce( Dm->Comm, nd.M); gnd.Px=sumReduce( Dm->Comm, nd.Px); gnd.Py=sumReduce( Dm->Comm, nd.Py); gnd.Pz=sumReduce( Dm->Comm, nd.Pz); gnd.K=sumReduce( Dm->Comm, nd.K); - gnd.p=sumReduce( Dm->Comm, nd.p); gwd.M=sumReduce( Dm->Comm, wd.M); gwd.Px=sumReduce( Dm->Comm, wd.Px); gwd.Py=sumReduce( Dm->Comm, wd.Py); gwd.Pz=sumReduce( Dm->Comm, wd.Pz); gwd.K=sumReduce( Dm->Comm, wd.K); - gwd.p=sumReduce( Dm->Comm, wd.p); gnc.M=sumReduce( Dm->Comm, nc.M); gnc.Px=sumReduce( Dm->Comm, nc.Px); gnc.Py=sumReduce( Dm->Comm, nc.Py); gnc.Pz=sumReduce( Dm->Comm, nc.Pz); gnc.K=sumReduce( Dm->Comm, nc.K); - gnc.p=sumReduce( Dm->Comm, nc.p); gwc.M=sumReduce( Dm->Comm, wc.M); gwc.Px=sumReduce( Dm->Comm, wc.Px); gwc.Py=sumReduce( Dm->Comm, wc.Py); gwc.Pz=sumReduce( Dm->Comm, wc.Pz); gwc.K=sumReduce( Dm->Comm, wc.K); - gwc.p=sumReduce( Dm->Comm, wc.p); giwn.Mn=sumReduce( Dm->Comm, iwn.Mn); giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx); @@ -624,4 +667,102 @@ void SubPhase::Full(){ } +void SubPhase::AggregateLabels(char *FILENAME){ + + int nx = Dm->Nx; + int ny = Dm->Ny; + int nz = Dm->Nz; + + int npx = Dm->nprocx(); + int npy = Dm->nprocy(); + int npz = Dm->nprocz(); + + int ipx = Dm->iproc(); + int ipy = Dm->jproc(); + int ipz = Dm->kproc(); + + int nprocs = Dm->nprocx()*Dm->nprocy()*Dm->nprocz(); + + int full_nx = npx*(nx-2); + int full_ny = npy*(ny-2); + int full_nz = npz*(nz-2); + int local_size = (nx-2)*(ny-2)*(nz-2); + long int full_size = long(full_nx)*long(full_ny)*long(full_nz); + + signed char *LocalID; + LocalID = new signed char [local_size]; + + //printf("aggregate labels: local size=%i, global size = %i",local_size, full_size); + // assign the ID for the local sub-region + for (int k=1; kid[n]; + if (local_id_val > 0){ + double value = Phi(i,j,k); + if (value > 0.0) local_id_val = 1; + else local_id_val = 2; + } + LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val; + } + } + } + MPI_Barrier(Dm->Comm); + + // populate the FullID + if (Dm->rank() == 0){ + signed char *FullID; + FullID = new signed char [full_size]; + // first handle local ID for rank 0 + for (int k=1; kComm,MPI_STATUS_IGNORE); + for (int k=1; krank(); + int dstrank = 0; + MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Dm->Comm); + } + MPI_Barrier(Dm->Comm); + +} + + diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h index 738d2188..683fc46a 100644 --- a/analysis/SubPhase.h +++ b/analysis/SubPhase.h @@ -101,7 +101,8 @@ public: void Basic(); void Full(); void Write(int time); - + void AggregateLabels(char *FILENAME); + private: FILE *TIMELOG; FILE *SUBPHASE; diff --git a/analysis/dcel.cpp b/analysis/dcel.cpp index 8267f09f..4c7be292 100644 --- a/analysis/dcel.cpp +++ b/analysis/dcel.cpp @@ -312,6 +312,7 @@ double DECL::EdgeAngle(int edge) } if (dotprod > 1.f) dotprod=1.f; + if (dotprod < -1.f) dotprod=-1.f; angle = acos(dotprod); /* project onto plane of cube face also works W = U - dotprod*V; @@ -344,6 +345,7 @@ double DECL::EdgeAngle(int edge) // concave angle = -angle; } + if (angle != angle) angle = 0.0; //printf("angle=%f,dot=%f (Edge=%i, twin=%i): P={%f, %f, %f}, Q={%f, %f, %f} U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),P.x,P.y,P.z,Q.x,Q.y,Q.z,U.x,U.y,U.z,V.x,V.y,V.z); return angle; } diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 4697aa20..8f25fad3 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -247,6 +247,9 @@ public: fillData.copy(Averages.Vel_z,VelzData); fillData.copy(Averages.morph_n->label,BlobData); IO::writeData( timestep, visData, comm.comm ); + char CurrentIDFilename[40]; + sprintf(CurrentIDFilename,"id_t%d.raw",timestep); + Averages.AggregateLabels(CurrentIDFilename); PROFILE_STOP("Save Vis",1); }; private: @@ -1020,4 +1023,3 @@ void runAnalysis::WriteVisData( int timestep, SubPhase &Averages, const double * PROFILE_STOP("write vis"); } - diff --git a/common/Domain.cpp b/common/Domain.cpp index 9ed0d89e..9c11d02a 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -128,6 +128,59 @@ Domain::Domain( std::shared_ptr db, MPI_Comm Communicator): rank_info = RankInfoStruct( myrank, rank_info.nx, rank_info.ny, rank_info.nz ); MPI_Barrier(Comm); } + +Domain::~Domain() +{ + // Free sendList + delete [] sendList_x; delete [] sendList_y; delete [] sendList_z; + delete [] sendList_X; delete [] sendList_Y; delete [] sendList_Z; + delete [] sendList_xy; delete [] sendList_yz; delete [] sendList_xz; + delete [] sendList_Xy; delete [] sendList_Yz; delete [] sendList_xZ; + delete [] sendList_xY; delete [] sendList_yZ; delete [] sendList_Xz; + delete [] sendList_XY; delete [] sendList_YZ; delete [] sendList_XZ; + // Free sendBuf + delete [] sendBuf_x; delete [] sendBuf_y; delete [] sendBuf_z; + delete [] sendBuf_X; delete [] sendBuf_Y; delete [] sendBuf_Z; + delete [] sendBuf_xy; delete [] sendBuf_yz; delete [] sendBuf_xz; + delete [] sendBuf_Xy; delete [] sendBuf_Yz; delete [] sendBuf_xZ; + delete [] sendBuf_xY; delete [] sendBuf_yZ; delete [] sendBuf_Xz; + delete [] sendBuf_XY; delete [] sendBuf_YZ; delete [] sendBuf_XZ; + // Free recvList + delete [] recvList_x; delete [] recvList_y; delete [] recvList_z; + delete [] recvList_X; delete [] recvList_Y; delete [] recvList_Z; + delete [] recvList_xy; delete [] recvList_yz; delete [] recvList_xz; + delete [] recvList_Xy; delete [] recvList_Yz; delete [] recvList_xZ; + delete [] recvList_xY; delete [] recvList_yZ; delete [] recvList_Xz; + delete [] recvList_XY; delete [] recvList_YZ; delete [] recvList_XZ; + // Free recvBuf + delete [] recvBuf_x; delete [] recvBuf_y; delete [] recvBuf_z; + delete [] recvBuf_X; delete [] recvBuf_Y; delete [] recvBuf_Z; + delete [] recvBuf_xy; delete [] recvBuf_yz; delete [] recvBuf_xz; + delete [] recvBuf_Xy; delete [] recvBuf_Yz; delete [] recvBuf_xZ; + delete [] recvBuf_xY; delete [] recvBuf_yZ; delete [] recvBuf_Xz; + delete [] recvBuf_XY; delete [] recvBuf_YZ; delete [] recvBuf_XZ; + // Free sendData + delete [] sendData_x; delete [] sendData_y; delete [] sendData_z; + delete [] sendData_X; delete [] sendData_Y; delete [] sendData_Z; + delete [] sendData_xy; delete [] sendData_xY; delete [] sendData_Xy; + delete [] sendData_XY; delete [] sendData_xz; delete [] sendData_xZ; + delete [] sendData_Xz; delete [] sendData_XZ; delete [] sendData_yz; + delete [] sendData_yZ; delete [] sendData_Yz; delete [] sendData_YZ; + // Free recvData + delete [] recvData_x; delete [] recvData_y; delete [] recvData_z; + delete [] recvData_X; delete [] recvData_Y; delete [] recvData_Z; + delete [] recvData_xy; delete [] recvData_xY; delete [] recvData_Xy; + delete [] recvData_XY; delete [] recvData_xz; delete [] recvData_xZ; + delete [] recvData_Xz; delete [] recvData_XZ; delete [] recvData_yz; + delete [] recvData_yZ; delete [] recvData_Yz; delete [] recvData_YZ; + // Free id + delete [] id; + // Free the communicator + if ( Comm != MPI_COMM_WORLD && Comm != MPI_COMM_NULL ) { + MPI_Comm_free(&Comm); + } +} + void Domain::initialize( std::shared_ptr db ) { d_db = db; @@ -187,58 +240,337 @@ void Domain::initialize( std::shared_ptr db ) MPI_Comm_size( Comm, &nprocs ); INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!"); } -Domain::~Domain() -{ - // Free sendList - delete [] sendList_x; delete [] sendList_y; delete [] sendList_z; - delete [] sendList_X; delete [] sendList_Y; delete [] sendList_Z; - delete [] sendList_xy; delete [] sendList_yz; delete [] sendList_xz; - delete [] sendList_Xy; delete [] sendList_Yz; delete [] sendList_xZ; - delete [] sendList_xY; delete [] sendList_yZ; delete [] sendList_Xz; - delete [] sendList_XY; delete [] sendList_YZ; delete [] sendList_XZ; - // Free sendBuf - delete [] sendBuf_x; delete [] sendBuf_y; delete [] sendBuf_z; - delete [] sendBuf_X; delete [] sendBuf_Y; delete [] sendBuf_Z; - delete [] sendBuf_xy; delete [] sendBuf_yz; delete [] sendBuf_xz; - delete [] sendBuf_Xy; delete [] sendBuf_Yz; delete [] sendBuf_xZ; - delete [] sendBuf_xY; delete [] sendBuf_yZ; delete [] sendBuf_Xz; - delete [] sendBuf_XY; delete [] sendBuf_YZ; delete [] sendBuf_XZ; - // Free recvList - delete [] recvList_x; delete [] recvList_y; delete [] recvList_z; - delete [] recvList_X; delete [] recvList_Y; delete [] recvList_Z; - delete [] recvList_xy; delete [] recvList_yz; delete [] recvList_xz; - delete [] recvList_Xy; delete [] recvList_Yz; delete [] recvList_xZ; - delete [] recvList_xY; delete [] recvList_yZ; delete [] recvList_Xz; - delete [] recvList_XY; delete [] recvList_YZ; delete [] recvList_XZ; - // Free recvBuf - delete [] recvBuf_x; delete [] recvBuf_y; delete [] recvBuf_z; - delete [] recvBuf_X; delete [] recvBuf_Y; delete [] recvBuf_Z; - delete [] recvBuf_xy; delete [] recvBuf_yz; delete [] recvBuf_xz; - delete [] recvBuf_Xy; delete [] recvBuf_Yz; delete [] recvBuf_xZ; - delete [] recvBuf_xY; delete [] recvBuf_yZ; delete [] recvBuf_Xz; - delete [] recvBuf_XY; delete [] recvBuf_YZ; delete [] recvBuf_XZ; - // Free sendData - delete [] sendData_x; delete [] sendData_y; delete [] sendData_z; - delete [] sendData_X; delete [] sendData_Y; delete [] sendData_Z; - delete [] sendData_xy; delete [] sendData_xY; delete [] sendData_Xy; - delete [] sendData_XY; delete [] sendData_xz; delete [] sendData_xZ; - delete [] sendData_Xz; delete [] sendData_XZ; delete [] sendData_yz; - delete [] sendData_yZ; delete [] sendData_Yz; delete [] sendData_YZ; - // Free recvData - delete [] recvData_x; delete [] recvData_y; delete [] recvData_z; - delete [] recvData_X; delete [] recvData_Y; delete [] recvData_Z; - delete [] recvData_xy; delete [] recvData_xY; delete [] recvData_Xy; - delete [] recvData_XY; delete [] recvData_xz; delete [] recvData_xZ; - delete [] recvData_Xz; delete [] recvData_XZ; delete [] recvData_yz; - delete [] recvData_yZ; delete [] recvData_Yz; delete [] recvData_YZ; - // Free id - delete [] id; - // Free the communicator - if ( Comm != MPI_COMM_WORLD && Comm != MPI_COMM_NULL ) { - MPI_Comm_free(&Comm); - } -} +void Domain::Decomp(std::shared_ptr domain_db ) +{ + //....................................................................... + // Reading the domain information file + //....................................................................... + int rank_offset = 0; + int RANK = rank(); + int nprocs, nprocx, nprocy, nprocz, nx, ny, nz; + int64_t global_Nx,global_Ny,global_Nz; + int64_t i,j,k,n; + int BC=0; + int64_t xStart,yStart,zStart; + int checkerSize; + //int inlet_layers_x, inlet_layers_y, inlet_layers_z; + //int outlet_layers_x, outlet_layers_y, outlet_layers_z; + xStart=yStart=zStart=0; + inlet_layers_x = 0; + inlet_layers_y = 0; + inlet_layers_z = 0; + outlet_layers_x = 0; + outlet_layers_y = 0; + outlet_layers_z = 0; + checkerSize = 32; + + // Read domain parameters + auto Filename = domain_db->getScalar( "Filename" ); + //auto L = domain_db->getVector( "L" ); + auto size = domain_db->getVector( "n" ); + auto SIZE = domain_db->getVector( "N" ); + auto nproc = domain_db->getVector( "nproc" ); + if (domain_db->keyExists( "offset" )){ + auto offset = domain_db->getVector( "offset" ); + xStart = offset[0]; + yStart = offset[1]; + zStart = offset[2]; + } + if (domain_db->keyExists( "InletLayers" )){ + auto InletCount = domain_db->getVector( "InletLayers" ); + inlet_layers_x = InletCount[0]; + inlet_layers_y = InletCount[1]; + inlet_layers_z = InletCount[2]; + } + if (domain_db->keyExists( "OutletLayers" )){ + auto OutletCount = domain_db->getVector( "OutletLayers" ); + outlet_layers_x = OutletCount[0]; + outlet_layers_y = OutletCount[1]; + outlet_layers_z = OutletCount[2]; + } + if (domain_db->keyExists( "checkerSize" )){ + checkerSize = domain_db->getScalar( "checkerSize" ); + } + else { + checkerSize = SIZE[0]; + } + auto ReadValues = domain_db->getVector( "ReadValues" ); + auto WriteValues = domain_db->getVector( "WriteValues" ); + auto ReadType = domain_db->getScalar( "ReadType" ); + + if (ReadType == "8bit"){ + } + else if (ReadType == "16bit"){ + } + else{ + //printf("INPUT ERROR: Valid ReadType are 8bit, 16bit \n"); + ReadType = "8bit"; + } + + nx = size[0]; + ny = size[1]; + nz = size[2]; + nprocx = nproc[0]; + nprocy = nproc[1]; + nprocz = nproc[2]; + global_Nx = SIZE[0]; + global_Ny = SIZE[1]; + global_Nz = SIZE[2]; + nprocs=nprocx*nprocy*nprocz; + char *SegData = NULL; + + if (RANK==0){ + + printf("Input media: %s\n",Filename.c_str()); + printf("Relabeling %lu values\n",ReadValues.size()); + for (int idx=0; idx 0){ + // use checkerboard pattern + printf("Checkerboard pattern at x inlet for %i layers \n",inlet_layers_x); + for (int k = 0; k 0){ + printf("Checkerboard pattern at y inlet for %i layers \n",inlet_layers_y); + // use checkerboard pattern + for (int k = 0; k 0){ + printf("Checkerboard pattern at z inlet for %i layers \n",inlet_layers_z); + // use checkerboard pattern + for (int k = zStart; k < zStart+inlet_layers_z; k++){ + for (int j = 0; j 0){ + // use checkerboard pattern + printf("Checkerboard pattern at x outlet for %i layers \n",outlet_layers_x); + for (int k = 0; k 0){ + printf("Checkerboard pattern at y outlet for %i layers \n",outlet_layers_y); + // use checkerboard pattern + for (int k = 0; k 0){ + printf("Checkerboard pattern at z outlet for %i layers \n",outlet_layers_z); + // use checkerboard pattern + for (int k = zStart + nz*nprocz - outlet_layers_z; k < zStart + nz*nprocz; k++){ + for (int j = 0; j LabelCount(ReadValues.size(),0); + // Set up the sub-domains + if (RANK==0){ + printf("Distributing subdomains across %i processors \n",nprocs); + printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz); + printf("Subdomain size: %i x %i x %i \n",nx,ny,nz); + printf("Size of transition region: %ld \n", z_transition_size); + + for (int kp=0; kp domain_db); void CommunicateMeshHalo(DoubleArray &Mesh); void CommInit(); int PoreCount(); diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 1654e81c..5e4f8eab 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -75,7 +75,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){ // read the input database db = std::make_shared( filename ); domain_db = db->getDatabase( "Domain" ); - color_db = db->getDatabase( "Color" ); + color_db = db->getDatabase( "Color" ); analysis_db = db->getDatabase( "Analysis" ); // set defaults @@ -166,8 +166,14 @@ void ScaLBL_ColorModel::SetDomain(){ } void ScaLBL_ColorModel::ReadInput(){ - size_t readID; - Mask->ReadIDs(); + + if (domain_db->keyExists( "Filename" )){ + Mask->Decomp(domain_db); + } + else{ + size_t readID; + Mask->ReadIDs(); + } for (int i=0; iid[i]; // save what was read sprintf(LocalRankString,"%05d",rank); @@ -648,6 +654,7 @@ void ScaLBL_ColorModel::Run(){ double volA = Averages->gnb.V; volA /= Dm->Volume; volB /= Dm->Volume;; + initial_volume = volA*Dm->Volume; double vA_x = Averages->gnb.Px/Averages->gnb.M; double vA_y = Averages->gnb.Py/Averages->gnb.M; double vA_z = Averages->gnb.Pz/Averages->gnb.M; @@ -681,10 +688,9 @@ void ScaLBL_ColorModel::Run(){ isSteady = true; if ( isSteady ){ - initial_volume = volA*Dm->Volume; MORPH_ADAPT = true; CURRENT_MORPH_TIMESTEPS=0; - delta_volume_target = Dm->Volume*morph_delta; //*volA ???? // set target volume change + delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change Averages->Full(); Averages->Write(timestep); analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); @@ -763,8 +769,7 @@ void ScaLBL_ColorModel::Run(){ else if (USE_SEED){ delta_volume = volA*Dm->Volume - initial_volume; CURRENT_MORPH_TIMESTEPS += analysis_interval; - double effective_seed_water = seed_water*(1.0+volB/volA); - double massChange = SeedPhaseField(effective_seed_water); + double massChange = SeedPhaseField(seed_water); if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target); } else if (USE_MORPHOPEN_OIL){ @@ -984,7 +989,6 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){ double mass_loss =0.f; double count =0.f; double *Aq_tmp, *Bq_tmp; - double SEED_THRESHOLD = -0.9; Aq_tmp = new double [7*Np]; Bq_tmp = new double [7*Np]; @@ -1012,58 +1016,54 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){ } } */ - double oil_value = 0.0; - double water_value = 1.0; for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){ + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; double phase_id = (dA - dB) / (dA + dB); - double random_value = double(rand())/ RAND_MAX; - if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){ - Aq_tmp[n] = 0.3333333333333333*oil_value; - Aq_tmp[n+Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value; + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; - Bq_tmp[n] = 0.3333333333333333*water_value; - Bq_tmp[n+Np] = 0.1111111111111111*water_value; - Bq_tmp[n+2*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+3*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+4*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+5*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+6*Np] = 0.1111111111111111*water_value; - mass_loss += oil_value - dA; - count++; + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; } + mass_loss += random_value*seed_water_in_oil; } for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){ - double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; - double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; + double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; double phase_id = (dA - dB) / (dA + dB); - double random_value = double(rand())/ RAND_MAX; - if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){ - Aq_tmp[n] = 0.3333333333333333*oil_value; - Aq_tmp[n+Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value; + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; - Bq_tmp[n] = 0.3333333333333333*water_value; - Bq_tmp[n+Np] = 0.1111111111111111*water_value; - Bq_tmp[n+2*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+3*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+4*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+5*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+6*Np] = 0.1111111111111111*water_value; - mass_loss += oil_value - dA; - count++; + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; } + mass_loss += random_value*seed_water_in_oil; } count= sumReduce( Dm->Comm, count); diff --git a/sample_scripts/configure_huckleberry b/sample_scripts/configure_huckleberry new file mode 100755 index 00000000..03d35824 --- /dev/null +++ b/sample_scripts/configure_huckleberry @@ -0,0 +1,30 @@ +# load the module for cmake +module load cmake +module load gcc/7.3.0 openmpi/3.1.2 cuda +module load szip +module load zlib +module list +module load hdf5/1.8.12 +module load silo + +cmake \ + -D CMAKE_C_COMPILER:PATH=mpicc \ + -D CMAKE_CXX_COMPILER:PATH=mpicxx \ + -D CMAKE_C_FLAGS="-fPIC" \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D MPI_COMPILER:BOOL=TRUE \ + -D MPIEXEC=mpirun \ + -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ + -D CMAKE_BUILD_TYPE:STRING=Release \ + -D USE_CUDA=1 \ + -D CMAKE_CUDA_FLAGS="-arch sm_60" \ + -D CMAKE_CUDA_HOST_COMPILER="$GCC_BIN/gcc" \ + -D USE_HDF5=1 \ + -D HDF5_DIRECTORY=$HDF5_DIR \ + -D USE_SILO=1 \ + -D SILO_DIRECTORY=$SILO_DIR \ + -D USE_TIMER=0 \ + ~/LBPM-WIA + +make VERBOSE=1 -j1 && make install + diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9737c7ad..b7362b09 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -39,6 +39,7 @@ CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_ # Add the tests ADD_LBPM_TEST( pmmc_cylinder ) +ADD_LBPM_TEST( TestSubphase ) ADD_LBPM_TEST( TestTorus ) ADD_LBPM_TEST( TestTorusEvolve ) ADD_LBPM_TEST( TestTopo3D ) @@ -67,6 +68,7 @@ ADD_LBPM_TEST_PARALLEL( TestSegDist 8 ) ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 ) ADD_LBPM_TEST_1_2_4( testCommunication ) ADD_LBPM_TEST( TestWriter ) +ADD_LBPM_TEST( TestDatabase ) IF ( USE_NETCDF ) ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 ) ADD_LBPM_EXECUTABLE( lbpm_uCT_pp ) diff --git a/tests/TestDatabase.cpp b/tests/TestDatabase.cpp new file mode 100644 index 00000000..58a35ee5 --- /dev/null +++ b/tests/TestDatabase.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include "common/UnitTest.h" +#include "common/Utilities.h" +#include "common/MPI_Helpers.h" +#include "common/Database.h" +#include "ProfilerApp.h" + + +// Main +int main(int argc, char **argv) +{ + int rank,nprocs; + MPI_Init(&argc,&argv); + MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); + Utilities::setAbortBehavior(true,2); + Utilities::setErrorHandlers(); + UnitTest ut; + int err=0; + + int BC=2; + int npx=1; int npy=2; int npz=4; + int nx=32; int ny=34; int nz=35; + std::vector List; + List.push_back("name1"); + List.push_back("name2"); + + // write a simple database and test that it can be read by LBPM + auto db = std::make_shared( ); + db->putScalar( "BC", BC ); + db->putVector( "nproc", { npx, npy, npz } ); + db->putVector( "n", { nx, ny, nz } ); + db->putVector( "Files", List); + + std::ofstream OutStream("test.db"); +// db->putDatabase(); + OutStream << "Domain { \n"; + db->print(OutStream, " "); + OutStream << "} \n"; + printf("TestDatbase: writing test file\n"); + OutStream.close(); + + auto new_db = std::make_shared( "test.db" ); + auto domain_db = new_db->getDatabase( "Domain" ); + if (domain_db->keyExists( "BC" )){ + auto newBC = domain_db->getScalar( "BC" ); + if (newBC != BC){ + err=1; + printf("TestDatbase: getScalar failed! \n"); + } + } + if (err==0) printf("TestDatabase: succeeded!\n"); + else printf("TestDatabase: errors encountered \n"); + + // Finished + PROFILE_SAVE("TestDatabase",true); + MPI_Barrier(comm); + MPI_Finalize(); + return err; +} + + diff --git a/tests/TestSubphase.cpp b/tests/TestSubphase.cpp new file mode 100644 index 00000000..8eb479bc --- /dev/null +++ b/tests/TestSubphase.cpp @@ -0,0 +1,146 @@ +// 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/SubPhase.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 for torus (Euler-Poincarie characteristic) \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)); + // const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + std::shared_ptr Averages(new SubPhase(Dm)); + Nx += 2; + Ny += 2; + Nz += 2; + int N = Nx*Ny*Nz; + //....................................................................... + for ( k=1;kid[n] = 1; + } + } + } + //....................................................................... + Dm->CommInit(); // Initialize communications for domains + //....................................................................... + + //....................................................................... + // Assign the phase ID field based and the signed distance + //....................................................................... + double R1,R2; + double CX,CY,CZ; //CY1,CY2; + CX=Nx*nprocx*0.5; + CY=Ny*nprocy*0.5; + CZ=Nz*nprocz*0.5; + R1 = Nx*nprocx*0.2; // middle radius + R2 = Nx*nprocx*0.1; // donut thickness + // + //CY1=Nx*nprocx*0.5+R1; + //CY2=Ny*nprocy*0.5-R1; + + double x,y,z; + if (rank==0) printf("Initializing the system \n"); + for ( k=1;kiproc()*Nx+i - CX; + y = Dm->jproc()*Ny+j - CY; + z = Dm->kproc()*Nz+k - CZ; + + // Shrink the sphere sizes by two voxels to make sure they don't touch + Averages->SDs(i,j,k) = 100.0; + //.............................................................................. + // Single torus + Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + // Double torus + /* y = Dm->jproc()*Ny+j - CY1; + //z = Dm->kproc()*Nz+k - CZ +R1; + Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2; + + y = Dm->jproc()*Ny+j - CY2; + //z = Dm->kproc()*Nz+k - CZ-R1; + Averages->Phi(i,j,k) = min(Averages->Phi(i,j,k), + sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2); + *///.............................................................................. + + //Averages->Phi(i,j,k) = - Averages->Phi(i,j,k); + if (Averages->Phi(i,j,k) > 0.0){ + Dm->id[n] = 2; + } + else{ + Dm->id[n] = 1; + } + } + } + } + + if (rank==0) printf("Aggregating labels \n"); + + Averages->AggregateLabels("phase_id.raw"); + // Averages->Reduce(); + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + MPI_Barrier(comm); + MPI_Finalize(); + return 0; +} + diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp index f73c63c5..a4b8e385 100644 --- a/tests/lbpm_serial_decomp.cpp +++ b/tests/lbpm_serial_decomp.cpp @@ -55,12 +55,16 @@ int main(int argc, char **argv) int64_t xStart,yStart,zStart; int checkerSize; int inlet_count_x, inlet_count_y, inlet_count_z; + int outlet_count_x, outlet_count_y, outlet_count_z; // char fluidValue,solidValue; xStart=yStart=zStart=0; inlet_count_x = 0; inlet_count_y = 0; inlet_count_z = 0; + outlet_count_x = 0; + outlet_count_y = 0; + outlet_count_z = 0; checkerSize = 32; // read the input database auto db = std::make_shared( filename ); @@ -84,6 +88,12 @@ int main(int argc, char **argv) inlet_count_y = InletCount[1]; inlet_count_z = InletCount[2]; } + if (domain_db->keyExists( "OutletLayers" )){ + auto OutletCount = domain_db->getVector( "OutletLayers" ); + outlet_count_x = OutletCount[0]; + outlet_count_y = OutletCount[1]; + outlet_count_z = OutletCount[2]; + } if (domain_db->keyExists( "checkerSize" )){ checkerSize = domain_db->getScalar( "checkerSize" ); } @@ -210,6 +220,63 @@ int main(int argc, char **argv) } } } + + if (outlet_count_x > 0){ + // use checkerboard pattern + printf("Checkerboard pattern at x outlet for %i layers \n",outlet_count_x); + for (int k = 0; k 0){ + printf("Checkerboard pattern at y outlet for %i layers \n",outlet_count_y); + // use checkerboard pattern + for (int k = 0; k 0){ + printf("Checkerboard pattern at z outlet for %i layers \n",outlet_count_z); + // use checkerboard pattern + for (int k = zStart + nz*nprocz - outlet_count_z; k < zStart + nz*nprocz; k++){ + for (int j = 0; j