Skip to content

add initial constitutive relations api - #17

Open
rrsettgast wants to merge 23 commits into
mainfrom
feature/addConstitutive
Open

add initial constitutive relations api#17
rrsettgast wants to merge 23 commits into
mainfrom
feature/addConstitutive

Conversation

@rrsettgast

Copy link
Copy Markdown
Contributor

This PR adds constitutive relations for activity and ionic strength.

constexpr IonicStrength::Params testParams
{
// Species charge
{ 1.0, -1.0, 2.0 }

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

The sum of charges needs to be 0. Suggesting to change 1.0 to -1.0

@frankfeifan

Copy link
Copy Markdown
Collaborator

Bdot tests are expected to fail now. Need to update the expected solution.

@frankfeifan frankfeifan self-assigned this Aug 14, 2026
Comment thread src/reactions/massActions/MassActions.hpp Outdated
Comment thread src/constitutive/activity/DebyeHuckel.hpp

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

I've heavily modified and refactored this file to account for the activity model during mass action solve. The root cause is secondary concentrations and activity coefficients are coupled (i.e., ionic strength depends on the secondary concentrations). So we should solve them together, and therefore, concentrations, activity coefficients and activities are mutually consistent.

The specific changes include:

  • Refactored calculateLogSecondarySpeciesConcentration and calculateLogSecondarySpeciesConcentrationWrtLogC: an inner Newton solve on log(C_sec) returning concentrations, activity coefficients, activities, a convergence flag, and the exact d log(C_sec)/d log(C_prim).

  • Renamed the old fixed-activity-coefficient versions to ...NoActivityUpdate and ...WrtLogCNoActivityUpdate.

  • Removed the two calculateAggregatePrimaryConcentrationsWrtLogC overloads and the calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC, which reconstructed secondary concentrations internally under gamma == 1. The two remaining aggregate functions are a pure mole balance, T_i = C_prim,i + sum_j nu_ji C_sec,j, taking the solved secondary arrays as input.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

New tests added:

Self-consistent solve, checking mass action holds and a = C * gamma for all species:

  • test_calculateLogSecondarySpeciesConcentration_allEquilibrium_identity
  • test_calculateLogSecondarySpeciesConcentration_allEquilibrium_bdot
  • test_calculateLogSecondarySpeciesConcentration_mixedSystem_identity
  • test_calculateLogSecondarySpeciesConcentration_mixedSystem_bdot
  • test_calculateLogSecondarySpeciesConcentration_identityActivityModel — reproduces the reference
    concentrations and confirms unit activity coefficients.

Derivative d log(C_sec)/d log(C_prim) against a central difference of the solve itself:

  • test_calculateLogSecondarySpeciesConcentrationWrtLogC_allEquilibrium_identity
  • test_calculateLogSecondarySpeciesConcentrationWrtLogC_allEquilibrium_bdot
  • test_calculateLogSecondarySpeciesConcentrationWrtLogC_mixedSystem_bdot

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants