diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index c753be18..960294fb 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -2544,6 +2544,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for x face(2,8,10,12,14)................................ ScaLBL_D3Q19_Pack(2, dvcSendList_x, 0, sendCount_x, sendbuf_x, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x, rank_x, sendtag + 0); req2[0] = @@ -2552,6 +2554,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for X face(1,7,9,11,13)................................ ScaLBL_D3Q19_Pack(1, dvcSendList_X, 0, sendCount_X, sendbuf_X, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X, rank_X, sendtag + 1); req2[1] = @@ -2560,6 +2564,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for y face(4,8,9,16,18)................................. ScaLBL_D3Q19_Pack(4, dvcSendList_y, 0, sendCount_y, sendbuf_y, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y, rank_y, sendtag + 2); req2[2] = @@ -2568,6 +2574,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for Y face(3,7,10,15,17)................................. ScaLBL_D3Q19_Pack(3, dvcSendList_Y, 0, sendCount_Y, sendbuf_Y, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y, rank_Y, sendtag + 3); req2[3] = @@ -2576,6 +2584,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for z face(6,12,13,16,17)................................ ScaLBL_D3Q19_Pack(6, dvcSendList_z, 0, sendCount_z, sendbuf_z, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z, rank_z, sendtag + 4); req2[4] = @@ -2584,6 +2594,8 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component) { //...Packing for Z face(5,11,14,15,18)................................ ScaLBL_D3Q19_Pack(5, dvcSendList_Z, 0, sendCount_Z, sendbuf_Z, &Aq[Component * 7 * N], N); + + ScaLBL_DeviceBarrier(); req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z, rank_Z, sendtag + 5); req2[5] = @@ -2697,6 +2709,7 @@ void ScaLBL_Communicator::TriSendD3Q7AA(double *Aq, double *Bq, double *Cq) { //................................................................................... // Send all the distributions + ScaLBL_DeviceBarrier(); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, 3 * sendCount_x, rank_x, sendtag + 0); req2[0] = @@ -2831,10 +2844,11 @@ void ScaLBL_Communicator::SendHalo(double *data) { ScaLBL_Scalar_Pack(dvcSendList_yZ, sendCount_yZ, sendbuf_yZ, data, N); ScaLBL_Scalar_Pack(dvcSendList_Yz, sendCount_Yz, sendbuf_Yz, data, N); ScaLBL_Scalar_Pack(dvcSendList_YZ, sendCount_YZ, sendbuf_YZ, data, N); + //................................................................................... // Send / Recv all the phase indcator field values //................................................................................... - + ScaLBL_DeviceBarrier(); req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x, rank_x, sendtag + 0); req2[0] = diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..282f5a11 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -111,6 +111,9 @@ void ScaLBL_ColorModel::ReadParams(string filename) { if (color_db->keyExists("flux")) { flux = color_db->getScalar("flux"); } + if (color_db->keyExists("timestep")) { + timestep = color_db->getScalar("timestep"); + } inletA = 1.f; inletB = 0.f; outletA = 0.f; @@ -516,7 +519,7 @@ void ScaLBL_ColorModel::Initialize() { double capillary_number = color_db->getScalar("capillary_number"); if (rank == 0) - printf(" set flux to achieve Ca=%f \n", capillary_number); + printf(" set flux to achieve Ca=%0.3e\n", capillary_number); double MuB = rhoB * (tauB - 0.5) / 3.0; double IFT = 6.0 * alpha; double CrossSectionalArea = @@ -524,7 +527,7 @@ void ScaLBL_ColorModel::Initialize() { flux = Mask->Porosity() * CrossSectionalArea * IFT * capillary_number / MuB; if (rank == 0) - printf(" flux=%f \n", flux); + printf(" flux=%0.3e\n", flux); } color_db->putScalar("flux", flux); @@ -620,8 +623,102 @@ void ScaLBL_ColorModel::Initialize() { ScaLBL_CopyToHost(Averages->Phi.data(), Phi, N * sizeof(double)); } +void ScaLBL_ColorModel::ForwardStep() { + // *************ODD TIMESTEP************* + timestep++; + // Compute the Phase indicator field + // Read for Aq, Bq happens in this routine (requires communication) + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + // Halo exchange for phase field + ScaLBL_Comm_Regular->SendHalo(Phi); + + ScaLBL_D3Q19_AAodd_Color( + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set BCs + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, + Fx, Fy, Fz, Nx, Nx * Ny, 0, + ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + + // *************EVEN TIMESTEP************* + timestep++; + // Compute the Phase indicator field + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); + + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + // Halo exchange for phase field + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + ScaLBL_Comm_Regular->SendHalo(Phi); + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, + Nx * Ny, ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, + Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + //************************************************************************ +} + double ScaLBL_ColorModel::Run(int returntime) { - int nprocs = nprocx * nprocy * nprocz; const RankInfoStruct rank_info(rank, nprocx, nprocy, nprocz); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); @@ -630,7 +727,6 @@ double ScaLBL_ColorModel::Run(int returntime) { bool RESCALE_FORCE = false; bool SET_CAPILLARY_NUMBER = false; bool TRIGGER_FORCE_RESCALE = false; - double tolerance = 0.01; auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); auto current_db = db->cloneDatabase(); auto flow_db = db->getDatabase("FlowAdaptor"); @@ -645,6 +741,12 @@ double ScaLBL_ColorModel::Run(int returntime) { double Ca_previous = 0.0; double minCa = 8.0e-6; double maxCa = 1.0; + + double force_mag_bracket_previous = 0.0; + double Ca_bracket_previous = 0.0; + double force_bracket_low, force_bracket_high, Ca_bracket_low; + bool bracket_found = false; + if (color_db->keyExists("capillary_number")) { capillary_number = color_db->getScalar("capillary_number"); SET_CAPILLARY_NUMBER = true; @@ -656,9 +758,8 @@ double ScaLBL_ColorModel::Run(int returntime) { color_db->getScalar("rescale_force_after_timestep"); RESCALE_FORCE = true; } - if (analysis_db->keyExists("tolerance")) { - tolerance = analysis_db->getScalar("tolerance"); - } + double tolerance = analysis_db->getWithDefault("tolerance", 0.01); + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map); @@ -666,105 +767,15 @@ double ScaLBL_ColorModel::Run(int returntime) { int CURRENT_TIMESTEP = 0; int EXIT_TIMESTEP = min(timestepMax, returntime); while (timestep < EXIT_TIMESTEP) { - PROFILE_START("Update"); - // *************ODD TIMESTEP************* - timestep++; - // Compute the Phase indicator field - // Read for Aq, Bq happens in this routine (requires communication) - ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, - ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, - ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - if (BoundaryCondition > 0 && BoundaryCondition < 5) { - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - // Halo exchange for phase field - ScaLBL_Comm_Regular->SendHalo(Phi); - - ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, - tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, - ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set BCs - if (BoundaryCondition == 3) { - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - if (BoundaryCondition == 4) { - din = - ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 5) { - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, - Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, - Fx, Fy, Fz, Nx, Nx * Ny, 0, - ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); + PROFILE_START("Update"); + ForwardStep(); - // *************EVEN TIMESTEP************* - timestep++; - // Compute the Phase indicator field - ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, - ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, - ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL - // Halo exchange for phase field - if (BoundaryCondition > 0 && BoundaryCondition < 5) { - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, - rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, - Nx * Ny, ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set boundary conditions - if (BoundaryCondition == 3) { - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 4) { - din = - ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 5) { - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, - rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, - Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); - //************************************************************************ analysis.basic( timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den); // allow initial ramp-up to get closer to steady state CURRENT_TIMESTEP += 2; - if (CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { + if (CURRENT_TIMESTEP % analysis_interval == 0 && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { analysis.finish(); double volB = Averages->gwb.V; @@ -801,7 +812,7 @@ double ScaLBL_ColorModel::Run(int returntime) { fabs(muA * flow_rate_A + muB * flow_rate_B) / (5.796 * alpha); bool isSteady = false; - if ((fabs((Ca - Ca_previous) / Ca) < tolerance && + if ((fabs((Ca - Ca_previous) / analysis_interval / Ca) < tolerance && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS)) isSteady = true; if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS) @@ -817,39 +828,72 @@ double ScaLBL_ColorModel::Run(int returntime) { timestep = INITIAL_TIMESTEP; TRIGGER_FORCE_RESCALE = true; if (rank == 0) - printf(" Capillary number missed target value = %f " - "(measured value was Ca = %f) \n ", + printf(" Capillary number missed target value = %0.3e " + "(measured value was Ca = %0.3e)\n", capillary_number, Ca); } - if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && - CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP) { + if (RESCALE_FORCE && SET_CAPILLARY_NUMBER && CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP) { TRIGGER_FORCE_RESCALE = true; } if (TRIGGER_FORCE_RESCALE) { RESCALE_FORCE = false; TRIGGER_FORCE_RESCALE = false; - double RESCALE_FORCE_FACTOR = capillary_number / Ca; - if (RESCALE_FORCE_FACTOR > 2.0) - RESCALE_FORCE_FACTOR = 2.0; - if (RESCALE_FORCE_FACTOR < 0.5) - RESCALE_FORCE_FACTOR = 0.5; - Fx *= RESCALE_FORCE_FACTOR; - Fy *= RESCALE_FORCE_FACTOR; - Fz *= RESCALE_FORCE_FACTOR; - force_mag = sqrt(Fx * Fx + Fy * Fy + Fz * Fz); - if (force_mag > 1e-3) { - Fx *= 1e-3 / force_mag; // impose ceiling for stability - Fy *= 1e-3 / force_mag; - Fz *= 1e-3 / force_mag; - } - if (rank == 0) - printf(" -- adjust force by factor %f \n ", - capillary_number / Ca); - Averages->SetParams(rhoA, rhoB, tauA, tauB, Fx, Fy, Fz, alpha, - beta); - color_db->putVector("F", {Fx, Fy, Fz}); + + if (bracket_found) { + // update brackets based on new results + if ((Ca - capillary_number) * (Ca_bracket_low - capillary_number) < 0) { + // on opposite side of target from low bracket, replace high bracket + force_bracket_high = force_mag; + } else { + // on same side of target as low bracket, replace low bracket + force_bracket_low = force_mag; + Ca_bracket_low = Ca; + } + } + + double RESCALE_FORCE_FACTOR; + if (!bracket_found) { + // if previous iteration was on opposite side of target Ca than current, then save bracket + if (force_mag_bracket_previous > 0 && (Ca_bracket_previous - capillary_number) * (Ca - capillary_number) < 0) { + force_bracket_low = std::min(force_mag_bracket_previous, force_mag); + force_bracket_high = std::max(force_mag_bracket_previous, force_mag); + Ca_bracket_low = (force_bracket_low == force_mag_bracket_previous) ? Ca_bracket_previous : Ca; + bracket_found = true; + } else + // still no bracket, fall back to proportional scaling + RESCALE_FORCE_FACTOR = std::min(std::max(capillary_number / Ca, 0.5), 2.0); + } + + if (bracket_found) + RESCALE_FORCE_FACTOR = 0.5 * (force_bracket_low + force_bracket_high) / force_mag; + + // save Ca and force for next iteration + force_mag_bracket_previous = force_mag; + Ca_bracket_previous = Ca; + + // rescale + if (rank == 0) + printf(" -- Setting rescale factor via %s: %0.3e\n", (bracket_found) ? "bracket" : "proportion", RESCALE_FORCE_FACTOR); + Fx *= RESCALE_FORCE_FACTOR; + Fy *= RESCALE_FORCE_FACTOR; + Fz *= RESCALE_FORCE_FACTOR; + force_mag = sqrt(Fx * Fx + Fy * Fy + Fz * Fz); + + if (force_mag > 1e-3 && Ca > minCa) { + if (rank == 0) + printf(" -- Body force ceiling reached. Target capillary number may be unreachable.\n"); + Fx *= 1e-3 / force_mag; // impose ceiling for stability + Fy *= 1e-3 / force_mag; + Fz *= 1e-3 / force_mag; + force_mag = 1e-3; + bracket_found = false; + } + if (rank == 0) + printf(" -- New force after rescaling, Fx: %0.3e, Fy: %0.3e, Fz: %0.3e\n", Fx, Fy, Fz); + Averages->SetParams(rhoA, rhoB, tauA, tauB, Fx, Fy, Fz, alpha, beta); + color_db->putVector("F", {Fx, Fy, Fz}); } if (isSteady) { Averages->Full(); @@ -859,8 +903,7 @@ double ScaLBL_ColorModel::Run(int returntime) { analysis.finish(); if (rank == 0) { - printf("** WRITE STEADY POINT *** "); - printf("Ca = %f, (previous = %f) \n", Ca, Ca_previous); + printf("Ca = %0.3e, (previous = %0.3e)\n", Ca, Ca_previous); double h = Dm->voxel_length; // pressures double pA = Averages->gnb.p; @@ -1048,7 +1091,7 @@ double ScaLBL_ColorModel::Run(int returntime) { } - printf(" Measured capillary number %f \n ", Ca); + printf(" Measured capillary number %0.3e\n", Ca); } if (SET_CAPILLARY_NUMBER) { Fx *= capillary_number / Ca; @@ -1060,18 +1103,27 @@ double ScaLBL_ColorModel::Run(int returntime) { Fz *= 1e-3 / force_mag; } if (rank == 0) - printf(" -- adjust force by factor %f \n ", + printf(" -- adjust force by factor %0.3e\n", capillary_number / Ca); Averages->SetParams(rhoA, rhoB, tauA, tauB, Fx, Fy, Fz, alpha, beta); color_db->putVector("F", {Fx, Fy, Fz}); + + // reset bisection search method + bracket_found = false; + force_mag_bracket_previous = 0.0; } else { if (rank == 0) { - printf("** Continue to simulate steady *** \n "); - printf("Ca = %f, (previous = %f) \n", Ca, Ca_previous); + printf("** Continue to simulate steady *** \n"); + printf("Ca = %0.3e, (previous = %0.3e)\n", Ca, Ca_previous); } } + + break; // steady-state achieved, exit. } + + // save for convergence checks + Ca_previous = Ca; } } analysis.finish(); @@ -1094,17 +1146,14 @@ double ScaLBL_ColorModel::Run(int returntime) { if (rank == 0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); return (MLUPS); - MLUPS *= nprocs; } void ScaLBL_ColorModel::Run() { int nprocs = nprocx * nprocy * nprocz; const RankInfoStruct rank_info(rank, nprocx, nprocy, nprocz); - int analysis_interval = - 1000; // number of timesteps in between in situ analysis - if (analysis_db->keyExists("analysis_interval")) { - analysis_interval = analysis_db->getScalar("analysis_interval"); - } + + int analysis_interval = analysis_db->getWithDefault("analysis_interval", 1000); + double tolerance = analysis_db->getWithDefault("tolerance", 0.0); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); @@ -1116,110 +1165,33 @@ void ScaLBL_ColorModel::Run() { Map); //analysis.createThreads( analysis_method, 4 ); auto t1 = std::chrono::system_clock::now(); - while (timestep < timestepMax) { - PROFILE_START("Update"); - // *************ODD TIMESTEP************* - timestep++; - // Compute the Phase indicator field - // Read for Aq, Bq happens in this routine (requires communication) - ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, - ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, - ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - if (BoundaryCondition > 0 && BoundaryCondition < 5) { - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - // Halo exchange for phase field - ScaLBL_Comm_Regular->SendHalo(Phi); - - ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, - tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, - ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set BCs - if (BoundaryCondition == 3) { - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - if (BoundaryCondition == 4) { - din = - ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 5) { - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, - Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, - Fx, Fy, Fz, Nx, Nx * Ny, 0, - ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); + double delta_sw = 1.0; + double sw_prev = -1.0; - // *************EVEN TIMESTEP************* - timestep++; - // Compute the Phase indicator field - ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, - ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, - ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL - // Halo exchange for phase field - if (BoundaryCondition > 0 && BoundaryCondition < 5) { - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, - rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, - Nx * Ny, ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set boundary conditions - if (BoundaryCondition == 3) { - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 4) { - din = - ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } else if (BoundaryCondition == 5) { - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, - rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, - Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); - //************************************************************************ + while (timestep < timestepMax && delta_sw > tolerance) { + PROFILE_START("Update"); + ForwardStep(); PROFILE_STOP("Update"); - if (rank == 0 && timestep % analysis_interval == 0 && - BoundaryCondition == 4) { - printf("%i %f \n", timestep, din); - } // Run the analysis analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den); + + if (timestep % analysis_interval == 0){ + analysis.finish(); + + double volA = Averages->gnb.V / Dm->Volume; + double volB = Averages->gwb.V / Dm->Volume; + double sw = volB / (volA + volB); + + delta_sw = fabs(sw - sw_prev) / analysis_interval / sw; + if (rank == 0) + printf("timestep: %d, sw: %0.5e, previous: %0.5e, delta: %.5e\n", + timestep, sw, sw_prev, delta_sw); + + sw_prev = sw; + } } analysis.finish(); PROFILE_STOP("Loop"); diff --git a/models/ColorModel.h b/models/ColorModel.h index 888f4fc2..975d984f 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -88,6 +88,11 @@ class ScaLBL_ColorModel { */ void Initialize(); + /** + * \brief Forward propagation step + */ + void ForwardStep(); + /** * \brief Run the simulation */ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ceafff4..ed5f3f0b 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -2,6 +2,7 @@ #ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator ) #ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator ) ADD_LBPM_EXECUTABLE( lbpm_color_simulator ) +ADD_LBPM_EXECUTABLE( lbpm_color_unsteady_simulator ) ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator ) diff --git a/tests/lbpm_color_unsteady_simulator.cpp b/tests/lbpm_color_unsteady_simulator.cpp new file mode 100644 index 00000000..67f91341 --- /dev/null +++ b/tests/lbpm_color_unsteady_simulator.cpp @@ -0,0 +1,65 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Utilities.h" +#include "models/ColorModel.h" + +/* + * Simulator for two-phase flow in porous media + * James E. McClure 2013-2014 + */ + + +//************************************************************************* +// Implementation of Two-Phase Immiscible LBM using CUDA +//************************************************************************* + +int main( int argc, char **argv ) +{ + + // Initialize + Utilities::startup( argc, argv, true ); + + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + + if ( rank == 0 ) { + printf( "********************************************************\n" ); + printf( "Running Color LBM Unsteady state \n" ); + printf( "********************************************************\n" ); + } + // Initialize compute device + int device = ScaLBL_SetDevice( rank ); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); + + PROFILE_ENABLE( 1 ); + PROFILE_SYNCHRONIZE(); + PROFILE_START( "Main" ); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); + ColorModel.ReadParams( filename ); + ColorModel.SetDomain(); + ColorModel.ReadInput(); + ColorModel.Create(); // creating the model will create data structure to match the pore + // structure and allocate variables + ColorModel.Initialize(); // initializing the model will set initial conditions for variables + + ColorModel.Run(); + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + cout << flush; + Utilities::shutdown(); + return 0; +}