diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 55c84560..790dd808 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -2186,101 +2186,291 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ } void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){ - - int Np = M.Np; - double dA, dB; - - double *Aq_tmp, *Bq_tmp; - - Aq_tmp = new double [7*Np]; - Bq_tmp = new double [7*Np]; - ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); - ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + int Np = M.Np; + double dA, dB; - for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ - 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]; - 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]; - if (dA > 1.0){ - double mass_change = dA - 1.0; - Aq_tmp[n] -= 0.333333333333333*mass_change; - Aq_tmp[n+Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - if (dB > 1.0){ - double mass_change = dB - 1.0; - Bq_tmp[n] -= 0.333333333333333*mass_change; - Bq_tmp[n+Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - } - for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ - 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]; - 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]; - if (dA > 1.0){ - double mass_change = dA - 1.0; - Aq_tmp[n] -= 0.333333333333333*mass_change; - Aq_tmp[n+Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - if (dB > 1.0){ - double mass_change = dB - 1.0; - Bq_tmp[n] -= 0.333333333333333*mass_change; - Bq_tmp[n+Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; - Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; - } - } + double *Aq_tmp, *Bq_tmp; - ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); - ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; + + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + + for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ + 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]; + 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]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ + 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]; + 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]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); } double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){ - - double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.1 ); - double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); - ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); - /* compute the local derivative of phase indicator field */ - double beta = M.beta; - double factor = 0.5/beta; - double total_interface_displacement = 0.0; - double total_interface_sites = 0.0; - for (int n=0; nPhi(n); - double dist1 = factor*log((1.0+value1)/(1.0-value1)); - double value2 = phi(n); - double dist2 = factor*log((1.0+value2)/(1.0-value2)); - phi_t(n) = value2; - if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ - /* time derivative of distance */ - double dxdt = 0.125*(dist2-dist1); - /* extrapolate to move the distance further */ - double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; - /* compute the new phase interface */ - phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); - total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); - total_interface_sites += 1.0; - } + double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.1 ); + double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); + + ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); + /* compute the local derivative of phase indicator field */ + double beta = M.beta; + double factor = 0.5/beta; + double total_interface_displacement = 0.0; + double total_interface_sites = 0.0; + for (int n=0; nPhi(n); + double dist1 = factor*log((1.0+value1)/(1.0-value1)); + double value2 = phi(n); + double dist2 = factor*log((1.0+value2)/(1.0-value2)); + phi_t(n) = value2; + if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ + /* time derivative of distance */ + double dxdt = 0.125*(dist2-dist1); + /* extrapolate to move the distance further */ + double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; + /* compute the new phase interface */ + phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); + total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); + total_interface_sites += 1.0; + } } - ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) ); - return total_interface_sites; + ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) ); + return total_interface_sites; +} + +double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){ + + const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz); + auto rank = M.rank; + auto Nx = M.Nx; auto Ny = M.Ny; auto Nz = M.Nz; + auto N = Nx*Ny*Nz; + double vF = 0.f; + double vS = 0.f; + double delta_volume; + double WallFactor = 1.0; + bool USE_CONNECTED_NWP = false; + + DoubleArray phase(Nx,Ny,Nz); + IntArray phase_label(Nx,Ny,Nz);; + DoubleArray phase_distance(Nx,Ny,Nz); + Array phase_id(Nx,Ny,Nz); + fillHalo fillDouble(M.Dm->Comm,M.Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1); + + // Basic algorithm to + // 1. Copy phase field to CPU + ScaLBL_CopyToHost(phase.data(), M.Phi, N*sizeof(double)); + + double count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f) count+=1.f; + } + } + } + double volume_initial = M.Dm->Comm.sumReduce( count); + double PoreVolume = M.Dm->Volume*M.Dm->Porosity(); + /*ensure target isn't an absurdly small fraction of pore volume */ + if (volume_initial < target_delta_volume*PoreVolume){ + volume_initial = target_delta_volume*PoreVolume; + } + + // 2. Identify connected components of phase field -> phase_label + + double volume_connected = 0.0; + double second_biggest = 0.0; + if (USE_CONNECTED_NWP){ + ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm); + M.Dm->Comm.barrier(); + + // only operate on component "0" + count = 0.0; + + for (int k=0; kComm.sumReduce( count); + second_biggest = M.Dm->Comm.sumReduce( second_biggest); + } + else { + // use the whole NWP + for (int k=0; kSDs(i,j,k) > 0.f){ + if (phase(i,j,k) > 0.f ){ + phase_id(i,j,k) = 0; + } + else { + phase_id(i,j,k) = 1; + } + } + else { + phase_id(i,j,k) = 1; + } + } + } + } + } + + // 3. Generate a distance map to the largest object -> phase_distance + CalcDist(phase_distance,phase_id,*M.Dm); + + double temp,value; + double factor=0.5/M.beta; + for (int k=0; k 1.f) value=1.f; + if (value < -1.f) value=-1.f; + // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. + temp = -factor*log((1.0+value)/(1.0-value)); + /// use this approximation close to the object + if (fabs(value) < 0.8 && M.Averages->SDs(i,j,k) > 1.f ){ + phase_distance(i,j,k) = temp; + } + // erase the original object + phase(i,j,k) = -1.0; + } + } + } + } + + + if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest ); + + if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial); + double target_delta_volume_incremental = target_delta_volume; + if (fabs(target_delta_volume) > 0.01*volume_initial) + target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume); + + delta_volume = MorphGrow(M.Averages->SDs,phase_distance,phase_id,M.Averages->Dm, target_delta_volume_incremental, WallFactor); + + for (int k=0; kSDs(i,j,k) > 0.f){ + if (d < 3.f){ + //phase(i,j,k) = -1.0; + phase(i,j,k) = (2.f*(exp(-2.f*M.beta*d))/(1.f+exp(-2.f*M.beta*d))-1.f); + } + } + } + } + } + fillDouble.fill(phase); + + count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f){ + count+=1.f; + } + } + } + } + double volume_final= M.Dm->Comm.sumReduce( count); + + delta_volume = (volume_final-volume_initial); + if (rank == 0) printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial); + if (rank == 0) printf(" new saturation = %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs))); + + // 6. copy back to the device + //if (rank==0) printf("MorphInit: copy data back to device\n"); + ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double)); + + // 7. Re-initialize phase field and density + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np); + auto BoundaryCondition = M.BoundaryCondition; + if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ + if (M.Dm->kproc()==0){ + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2); + } + if (M.Dm->kproc() == M.nprocz-1){ + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3); + } + } + return delta_volume; } diff --git a/models/ColorModel.h b/models/ColorModel.h index becb528b..2b76ef2f 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -99,6 +99,7 @@ public: ~FlowAdaptor(); double MoveInterface(ScaLBL_ColorModel &M); double ImageInit(ScaLBL_ColorModel &M, std::string Filename); + double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume); double UpdateFractionalFlow(ScaLBL_ColorModel &M); void Flatten(ScaLBL_ColorModel &M); DoubleArray phi; diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 3591bbcd..e67c48f8 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -86,7 +86,7 @@ int main( int argc, char **argv ) if (ColorModel.db->keyExists( "FlowAdaptor" )){ auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" ); MAX_STEADY_TIME = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); - SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 100000 ); + SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 50000 ); FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault( "fractional_flow_increment", 0.05); ENDPOINT_THRESHOLD = flow_db->getWithDefault( "endpoint_threshold", 0.1); } @@ -109,6 +109,7 @@ int main( int argc, char **argv ) /* update the fluid configuration with the flow adapter */ int skip_time = 0; timestep = ColorModel.timestep; + /* get the averaged flow measures computed internally for the last simulation point*/ double SaturationChange = 0.0; double volB = ColorModel.Averages->gwb.V; double volA = ColorModel.Averages->gnb.V; @@ -121,6 +122,7 @@ int main( int argc, char **argv ) double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); + /* stop simulation if previous point was sufficiently close to the endpoint*/ if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false; if (ContinueSimulation){ while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){ @@ -128,8 +130,11 @@ int main( int argc, char **argv ) if (PROTOCOL == "fractional flow") { Adapt.UpdateFractionalFlow(ColorModel); } + else if (PROTOCOL == "shell aggregation"){ + double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange; + Adapt.ShellAggregation(ColorModel,target_volume_change); + } else if (PROTOCOL == "image sequence"){ - // Use image sequence IMAGE_INDEX++; if (IMAGE_INDEX < IMAGE_COUNT){ std::string next_image = ImageList[IMAGE_INDEX];