Skip to content

bugfix: min_steady_timesteps in lbpm_color_simulator #107

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jeremyfirst22 wants to merge 1 commit into
base: dev
Choose a base branch
from
Open

bugfix: min_steady_timesteps in lbpm_color_simulator #107

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions models/ColorModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -630,7 +630,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<std::string>( "WettingConvention", "none" );
auto current_db = db->cloneDatabase();
auto flow_db = db->getDatabase("FlowAdaptor");
Expand All @@ -656,9 +655,8 @@ double ScaLBL_ColorModel::Run(int returntime) {
color_db->getScalar<int>("rescale_force_after_timestep");
RESCALE_FORCE = true;
}
if (analysis_db->keyExists("tolerance")) {
tolerance = analysis_db->getScalar<double>("tolerance");
}
double tolerance = analysis_db->getWithDefault<double>("tolerance", 1e-5);
int analysis_interval = analysis_db->getWithDefault<int>("analysis_interval", 1000);

runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular,
Map);
Expand Down Expand Up @@ -764,7 +762,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
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;
Expand Down Expand Up @@ -801,7 +799,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))
Copy link
Collaborator

@JamesEMcClure JamesEMcClure Nov 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why add analysis_interval here? The previous logic is simply comparing the relative change in the capillary number. This change will alter the interpretation of tolerance

Copy link
Contributor Author

@jeremyfirst22 jeremyfirst22 Dec 1, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JamesEMcClure I added the analysis_interval here to make the tolerance relative to each simulation step, and independent of the number of steps between evaluating the capillary number.

i.e., if the simulation is not changing by greater than 1% per step, we consider the simulation converged. Not if the simulation is changing by 1% per analysis interval.

The idea was that you should be free to change the analysis_interval without having to adjust the tolerance in kind to achieve the same results.

I agree this changes the interpretation of tolerance, but since the Ca_previous was never updated anyway (i.e., tolerance was never used), I figured I was free to update the interpretation. But 1% per step is an extremely loose default. If you'd like, I will update the default to be consistent with the previous interpretation (1% / 1000 steps)

Copy link
Contributor Author

@jeremyfirst22 jeremyfirst22 Dec 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JamesEMcClure I went ahead and make this change, and updated the default tolerance to 1e-5 so that is consistent with the old interpretation of the default (i..e, 0.01 over the default 1000 steps).

Let me know if you have any other comments on this one.

isSteady = true;
if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS)
Expand Down Expand Up @@ -1071,7 +1069,12 @@ double ScaLBL_ColorModel::Run(int returntime) {
printf("Ca = %f, (previous = %f) \n", Ca, Ca_previous);
}
}

break; // steady-state achieved, exit.
}

// save for convergence checks
Ca_previous = Ca;
}
}
analysis.finish();
Expand Down Expand Up @@ -1100,11 +1103,9 @@ double ScaLBL_ColorModel::Run(int returntime) {
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<int>("analysis_interval");
}

int analysis_interval = analysis_db->getWithDefault<int>("analysis_interval", 1000);
double tolerance = analysis_db->getWithDefault<double>("tolerance", 0.0);

//************ MAIN ITERATION LOOP ***************************************/
comm.barrier();
Expand Down