body force adjustment: scale-then-bisect algorithm

[jeremyfirst22]

Feature

• Improved body force search algorithm for steady state calculations
• Currently, the body force is proportionally adjusted to achieve the desired capillary number, but because these properties are non-linearly related, the proportional scaling was prone to 'straddling' the target and infinitely bouncing back-and-forth (as identified in Endless overshoot/undershoot Capillary number using Fractional Flow Protocol #75).
• Now, once the proportional scaling brackets the target, we switch to the bisection method for guaranteed convergence.
• Resolves Endless overshoot/undershoot Capillary number using Fractional Flow Protocol #75.

To do

• Investigating potential further issue for low capillary number simulations, where the desired capillary number is below that of the spurious current. This results in an 'infinite halving' of the body force, where each lowering does not affect the measured capillary number since it is already dominated by spurious currents.

[JamesEMcClure]

The main comment that I have with this is about how to test this feature (without having much background on how you may have tested it so far):

One of the main challenges with the algorithm to adjust the capillary number is that it can take a lot of timesteps for the flow velocity to reach a steady-state based on the driving force. If you adjust the force before the flow velocity has converged then you end up correcting based on a hasty prediction, in which case you will end up adjusting the force more often than necessary, which slows things down. The number of timesteps required to reach steady state will increase with system size. To confirm that this update works as intended I would run steady-state relative permeability simulations (e.g. you can use the random close sphere pack, e.g. example/Sph125 or example/Sph1896, or any representative uCT image) for a range of system sizes: 256^3, 512^3, 1024^3 .

[jeremyfirst22]

@JamesEMcClure I've tested this on the Sph125 as you suggested, but it seems the sphere pack capillary number response is too well-behaved to even trigger the bracketing approach; the capillary number target is reached just fine via the proportional scaling.

But on a number of real rocks, you get convergence behavior like this:

    Capillary number missed target value = 3.612e-06 (measured value was Ca = 4.273e-05)
    -- Setting rescale factor via proportion: 5.000e-01
    -- New force after rescaling, Fx: 0.000e+00, Fy: 0.000e+00, Fz: -5.000e-06
    Capillary number missed target value = 3.612e-06 (measured value was Ca = 1.984e-05)
    -- Setting rescale factor via proportion: 5.000e-01
    -- New force after rescaling, Fx: 0.000e+00, Fy: 0.000e+00, Fz: -2.500e-06
    Capillary number missed target value = 3.612e-06 (measured value was Ca = 8.061e-06)
    -- Setting rescale factor via proportion: 5.000e-01
    -- New force after rescaling, Fx: 0.000e+00, Fy: 0.000e+00, Fz: -1.250e-06
    Capillary number missed target value = 3.612e-06 (measured value was Ca = 2.536e-06)
    -- Setting rescale factor via bracket: 1.500e+00
    -- New force after rescaling, Fx: 0.000e+00, Fy: 0.000e+00, Fz: -1.875e-06
Ca = 5.464e-06, (previous = 2.536e-06)
  Measured capillary number 5.464e-06
    -- adjust force by factor 6.611e-01

Once the capillary number is overshot by the proportional scaling, it switches over to bisectional search as you can see in the output.

I've also added a switch to the convergence algorithm so that if the capillary number is no longer responding to the body force adjustment, print out a warning like this:

    WARNING: Target capillary number (3.612e-06) is unreachable. Continuing at current (9.129e-06)

and continue. I've seen this pop up in two cases, where the target capillary number is below spurious current floor (which is independent of this PR and you can test this in Sph125 by setting target capillary number to 1e-6), and when the capillary number response is quite noisy so the bracket collapses to a single body force, but still misses the target range.

There are likely more corner cases to catch and fix, but this has reduced the number of 'infinitely looping' cases I've seen pop on our cluster.