Reformulate solubility_products constraints

[lxhowl]

Fixes

pH-dependent solubility tests from WaterTAP

Summary/Motivation:

This solubility equilibrium constraint used a smoothing strategy to compare precipitate amount s and the solution state Q. Only one of s and Q can be greater than zero at any time. When s and Q differ in several orders of magnitude, the problem becomes hard to scale, especially when Q is in log form, the original formulation will sacrifice the accuracy of s, the precipitate amount, which could be very small (i.e., 1e-10) and important in water-related applications.

Fix outside of bounds warning in the initialization of pH-dependent solubility tests.

Changes proposed in this PR:

• Normalize the precipitate amount s to be 0-1 with respect to its solubility product constant Ksp.

• Update mole_frac_phase_comp bounds in FpcTP to be consistent with log_mole_frac_phase_comp bounds.

Legal Acknowledgement

By contributing to this software project, I agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the license terms described in the LICENSE.txt file at the top level of this directory.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[bknueven]

For a bit more context, Xinhong has been working to make some of chemistry models in WaterTAP more stable. She found this issue with the formulation and came up with this fix to keep these smooth_max constraints scaled consistently.

[andrewlee94]

I have a few questions about why you chose the formulation you did, as well as about consitency with other state definitions.

[andrewlee94]

Also, the tests for the solubility product are all failing right now.

[adowling2]

An alternate approach is to add a scaling factor to smooth_max.

scaled_smooth_max(0,x,scale) = scale*smooth_max(0, x/scale)

From an implementation perspective, I proposed adding an optional scale=1 argument to smooth_max that does this.

I also like @andrewlee94's suggestion about scaling variables carefully.

[andrewlee94]

@lxhowl Changing the bounds on the mole fractions is not the correct way to address the variables outside of bounds warning - that is something that needs to be fixed in either the model or the tool that is setting the value outside the bounds. The value of 1w-20 was deliberately chose to avoid singularities, and this change is likely what caused the testing in the plate heat exchanger model to start failing.

[codecov]

Codecov Report

Patch coverage: 84.61% and project coverage change: -0.01% ⚠️

Comparison is base (ad5b8ad) 76.84% compared to head (14cb740) 76.83%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1197      +/-   ##
==========================================
- Coverage   76.84%   76.83%   -0.01%     
==========================================
  Files         390      390              
  Lines       61875    61885      +10     
  Branches    11389    11391       +2     
==========================================
+ Hits        47547    47552       +5     
- Misses      11866    11871       +5     
  Partials     2462     2462              

Files Changed	Coverage Δ
.../modular_properties/reactions/equilibrium_forms.py	87.21% <84.61%> (-0.59%)	⬇️

... and 1 file with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[agarciadiego]

@lxhowl if the user does not specify k_eq_ref how robust is the s_norm default value? Have we tested this with another example where the k_eq_ref is unknown? @andrewlee94 would it make more sense to have an if statement to use this formulation in case the user has an a priori knowledge of k_eq_ref

[andrewlee94]

@agarciadiego I am leaning towards just having the one version - the old one already had a smoothing factor that the user should be setting anyway so this just adds more flexibility (and better scaling overall).

If we think it is necessary to keep both formulations, then I think they should be separate methods the user chooses from rather than an if statement in the code. An if buried in the code is far harder for the user to find and understand.

@lxhowl This discussion did remind me however that you need to update the documentation to reflect the new formulation and explain the new parameters (and how the user should adjust them).

[lxhowl]

@agarciadiego The tests I updated contains cases without k_eq_ref specified. I will add assertions for s_norm and s_scale default values checking.

The formulation of the updated s term is meant to keep it in the same level of accuracy as the Q term. It is doing a scaling down when s is large, so that we can preserve accuracy of Q term. And when s is small, it scale up to keep higher accuracy of s. As s_norm decreases, very low concentration s can be captured. The default s_norm value 1e-4 without k_eq_ref specified is set in considerations of common Ksp range 1e-37~1e-4.

[lbianchi-lbl]

The test failures are likely due to #1215

[eslickj]

If one or both Q and s need to be zero, could this just bee smooth min(Q, s) == 0?

[lxhowl]

Current expression goes Q-0.5(Q-s+|Q-s|) == 0 while the smooth min gives 0.5(Q+s-|Q-s|)==0. The complexity remains the same. The current expression keeps Q as a leading term, which might be easy for users to scale the solubility product with other equilibrium reactions.

[andrewlee94]

@eslickj I think this formulation is correct - it is a bit more complicated than just one of Q or s needing to be zero.

[eslickj]

@andrewlee94, I'm sure you are probably right. Can someone explain the difference to me? Why is this one preferred? I think both should work, but I may be missing something. To me, the the one in the code currently just seems convoluted.

[andrewlee94]

@eslickj I cannot remember the details - I asked Larry Biegler to help with this is and this is the form he suggested.

[andrewlee94]

@lxhowl Is this ready for re-review yet?

[lxhowl]

@lxhowl Is this ready for re-review yet?
YES, it is ready.

[andrewlee94]

Looks good to me.