GeochemistryActivityCoefficientsDebyeHuckelTest.C File Reference

Go to the source code of this file.

Functions
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, notDH)
	Test errors when a non-DebyeHuckel database is read. More...
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, noNeutral)
	Test behavior when a database with incorrect neutral species is read. More...
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, ionicStrengthInterface)
	Test interface with ionic strength object. More...
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, getDebyeHuckel)
	Test getDebyeHuckel. More...
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, buildActivityCoefficientsDebyeHuckel)
	Test calculate activity coefficients for method=DebyeHuckel. More...
	TEST (GeochemistryActivityCoefficientsDebyeHuckelTest, waterActivity)
	Test water activity for method=DebyeHuckel. More...

Variables
const GeochemicalDatabaseReader	database ("database/moose_testdb.json", true, true, false)
const PertinentGeochemicalSystem	model (database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
const ModelGeochemicalDatabase	mgd = model.modelGeochemicalDatabase()

Function Documentation

TEST() [1/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		notDH
	)

Test errors when a non-DebyeHuckel database is read.

Definition at line 31 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  const GeochemicalDatabaseReader database_notDH(
      "database/faultydbs/notDH.json", true, true, false);
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  try
  {
    GeochemistryActivityCoefficientsDebyeHuckel ac_notDH(is, database_notDH);
    FAIL() << "Missing expected exception.";
  }
  catch (const std::exception & e)
  {
    std::string msg(e.what());
    ASSERT_TRUE(msg.find("Attempted to get Debye-Huckel activity parameters but the activity model "
                         "is NOTdebye-huckel") != std::string::npos)
        << "Failed with unexpected error message: " << msg;
  }
}

TEST() [2/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		noNeutral
	)

Test behavior when a database with incorrect neutral species is read.

Definition at line 51 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  const GeochemicalDatabaseReader database_noNeutral(
      "database/faultydbs/missing_DH.json", true, true, false);
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  GeochemistryActivityCoefficientsDebyeHuckel ac_noNeutral(is, database_noNeutral);
  ac_noNeutral.setInternalParameters(
      25.0, mgd, std::vector<Real>(8, 1.0), std::vector<Real>(9), std::vector<Real>(3));
  const DebyeHuckelParameters dh = ac_noNeutral.getDebyeHuckel();

  EXPECT_EQ(dh.A, 0.5555);
  EXPECT_EQ(dh.B, 0.3333);
  EXPECT_EQ(dh.Bdot, 0.0444);
  EXPECT_EQ(dh.a_water, 0.0);
  EXPECT_EQ(dh.b_water, 0.0);
  EXPECT_EQ(dh.c_water, 0.0);
  EXPECT_EQ(dh.d_water, 0.0);
  EXPECT_EQ(dh.a_neutral, 0.0);
  EXPECT_EQ(dh.b_neutral, 0.0);
  EXPECT_EQ(dh.c_neutral, 0.0);
  EXPECT_EQ(dh.d_neutral, 0.0);
}

TEST() [3/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		ionicStrengthInterface
	)

Test interface with ionic strength object.

Definition at line 75 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  GeochemistryActivityCoefficientsDebyeHuckel ac(is, database);

  is.setUseOnlyBasisMolality(true);
  is.setMaxStoichiometricIonicStrength(2.0E-8);

  Real gold_ionic_str = 0.0;
  std::vector<Real> basis_m(8);
  basis_m[mgd.basis_species_index.at("H2O")] = 1.0;
  basis_m[mgd.basis_species_index.at("H+")] = 2.0;
  gold_ionic_str += 2.0;
  basis_m[mgd.basis_species_index.at("HCO3-")] = 3.0;
  gold_ionic_str += 3.0;
  basis_m[mgd.basis_species_index.at("O2(aq)")] = 4.0;
  basis_m[mgd.basis_species_index.at("Ca++")] = 5.0;
  gold_ionic_str += 4 * 5.0;
  basis_m[mgd.basis_species_index.at(">(s)FeOH")] = 6.0;
  basis_m[mgd.basis_species_index.at("radius_neg1")] = 7.0;
  basis_m[mgd.basis_species_index.at("radius_neg1.5")] = 8.0;

  // the equilibrium and kinetic molalities do not make any difference due to ionic strength using
  // only basis molalities
  ac.setInternalParameters(
      25.0, mgd, basis_m, std::vector<Real>(9, 1.1), std::vector<Real>(3, 2.2));

  EXPECT_NEAR(ac.getIonicStrength(), 0.5 * gold_ionic_str, 1E-8);
  EXPECT_EQ(ac.getStoichiometricIonicStrength(), 2.0E-8);
}

TEST() [4/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		getDebyeHuckel
	)

Test getDebyeHuckel.

Definition at line 107 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  GeochemistryActivityCoefficientsDebyeHuckel ac(is, database);
  ac.setInternalParameters(
      25.0, mgd, std::vector<Real>(8, 1.0), std::vector<Real>(9), std::vector<Real>(3));
  const DebyeHuckelParameters dh = ac.getDebyeHuckel();

  EXPECT_EQ(dh.A, 0.5092);
  EXPECT_EQ(dh.B, 0.3283);
  EXPECT_EQ(dh.Bdot, 0.041);
  EXPECT_EQ(dh.a_water, 1.45397);
  EXPECT_EQ(dh.b_water, 0.022357);
  EXPECT_EQ(dh.c_water, 0.0093804);
  EXPECT_EQ(dh.d_water, -0.0005362);
  EXPECT_EQ(dh.a_neutral, 0.1127);
  EXPECT_EQ(dh.b_neutral, -0.01049);
  EXPECT_EQ(dh.c_neutral, 0.001545);
  EXPECT_EQ(dh.d_neutral, 0.0);
}

TEST() [5/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		buildActivityCoefficientsDebyeHuckel
	)

Test calculate activity coefficients for method=DebyeHuckel.

Definition at line 129 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  GeochemistryActivityCoefficientsDebyeHuckel ac(is, database);
  ac.setInternalParameters(
      25.0, mgd, std::vector<Real>(8, 1.0), std::vector<Real>(9), std::vector<Real>(3));

  const Real ionic_str = ac.getIonicStrength();

  std::vector<Real> basis_ac;
  std::vector<Real> eqm_ac;
  ac.buildActivityCoefficients(mgd, basis_ac, eqm_ac);

  EXPECT_EQ(basis_ac.size(), (std::size_t)8);
  EXPECT_EQ(eqm_ac.size(), (std::size_t)9);

  const DebyeHuckelParameters dh = ac.getDebyeHuckel();

  // Note that GeochemistryActivity has been tested elsewhere: here we're just checking the slots
  // get filled correctly
  for (unsigned i = 1; i < 8; ++i) // don't loop over water
  {
    Real gold = 0.0;
    if (i == 3) // O2(aq)
      gold = std::pow(10.0,
                      GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
                          ionic_str, dh.a_neutral, dh.b_neutral, dh.c_neutral, dh.d_neutral));
    else if (i == 5) // >(s)FeOH
      gold = 1.0;
    else if (i == 6) // radius_neg1
      gold = std::pow(
          10.0,
          GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(ionic_str, dh.Bdot));
    else if (i == 7) // radius_neg1.5
      gold = 1.0;
    else
      gold =
          std::pow(10.0,
                   GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.basis_species_charge[i],
                                                                        mgd.basis_species_radius[i],
                                                                        std::sqrt(ionic_str),
                                                                        dh.A,
                                                                        dh.B,
                                                                        dh.Bdot));
    ASSERT_NEAR(basis_ac[i], gold, 1.0E-8);
  }

  for (unsigned j = 0; j < 9; ++j)
  {
    Real gold = 0.0;
    const std::string nm = mgd.eqm_species_name[j];
    if (nm == "Calcite")
      continue;
    else if (nm == "CO2(aq)" || nm == "CaCO3" || nm == ">(s)FeO-")
      gold = 1.0;
    else if (nm == "seq_radius_neg1")
      gold = std::pow(
          10.0,
          GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(ionic_str, dh.Bdot));
    else
      gold =
          std::pow(10.0,
                   GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.eqm_species_charge[j],
                                                                        mgd.eqm_species_radius[j],
                                                                        std::sqrt(ionic_str),
                                                                        dh.A,
                                                                        dh.B,
                                                                        dh.Bdot));
    ASSERT_NEAR(eqm_ac[j], gold, 1.0E-8);
  }
}

TEST() [6/6]

TEST	(	GeochemistryActivityCoefficientsDebyeHuckelTest	,
		waterActivity
	)

Test water activity for method=DebyeHuckel.

Definition at line 202 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

{
  GeochemistryIonicStrength is(1.0E9, 2.0E9, false, false);
  GeochemistryActivityCoefficientsDebyeHuckel ac(is, database);
  ac.setInternalParameters(
      25.0, mgd, std::vector<Real>(8, 1.0), std::vector<Real>(9), std::vector<Real>(3));

  const DebyeHuckelParameters dh = ac.getDebyeHuckel();
  const Real stoi_ionic_str = ac.getStoichiometricIonicStrength();

  // note that GeochemistryActivityCalculators::lnActivityDHBdotWater has been tested elsewhere:
  // this is just testing the information gets passed through to
  // GeochemistryActivityCoefficientsDebyeHuckel
  EXPECT_NEAR(ac.waterActivity(),
              std::exp(GeochemistryActivityCalculators::lnActivityDHBdotWater(
                  stoi_ionic_str, dh.A, dh.a_water, dh.b_water, dh.c_water, dh.d_water)),
              1.0E-8);
}

Variable Documentation

database

const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)

Referenced by TEST(), TEST_F(), and testExceptionMessage().

mgd

const ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase()

Definition at line 28 of file GeochemistryActivityCoefficientsDebyeHuckelTest.C.

Referenced by AddGeochemistrySolverAction::act(), GeochemistrySpeciesSwapper::alterMGD(), GeochemistryActivityCoefficientsDebyeHuckel::buildActivityCoefficients(), GeochemistrySpeciesSwapper::checkSwap(), GeochemistryQuantityAux::computeValue(), GeochemistrySpeciesSwapper::constructInverseMatrix(), GeochemistrySpeciesSwapper::findBestEqmSwap(), GeochemistryUnitConverter::fromMoles(), GeochemistryQuantityAux::GeochemistryQuantityAux(), GeochemistryIonicStrength::ionicStrength(), GeochemistryConsoleOutput::output(), GeochemistrySpeciesSwapper::performSwap(), GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters(), GeochemicalSystem::setModelGeochemicalDatabase(), GeochemicalSolver::solveSystem(), GeochemistryIonicStrength::stoichiometricIonicStrength(), GeochemicalSolver::swapNeeded(), TEST(), TEST_F(), and GeochemistryUnitConverter::toMoles().

model

const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")

Referenced by MultiPlasticityRawComponentAssembler::activeModelSurfaces(), MultiPlasticityRawComponentAssembler::activeSurfaces(), MultiPlasticityRawComponentAssembler::anyActiveSurfaces(), MultiPlasticityRawComponentAssembler::buildActiveConstraints(), MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), MultiPlasticityDebugger::checkSolution(), ModularGapConductanceConstraint::computeQpResidual(), ADComputeSmearedCrackingStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityDebugger::fddflowPotential_dintnl(), MultiPlasticityDebugger::fddyieldFunction_dintnl(), MultiPlasticityDebugger::fdJacobian(), MultiPlasticityRawComponentAssembler::flowPotential(), ExtremeValue::function(), Gaussian::function(), TruncatedGaussian::function(), MultiPlasticityRawComponentAssembler::hardPotential(), CombinedScalarDamageTempl< is_ad >::initialSetup(), ComputeMultipleCrystalPlasticityStress::initialSetup(), LoadCovarianceDataAction::load(), MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), MultiPlasticityLinearSystem::nrStep(), MultiPlasticityDebugger::outputAndCheckDebugParameters(), ComputeMultiPlasticityStress::plasticStep(), ComputeMultiPlasticityStress::residual2(), ComputeMultiPlasticityStress::returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), ComputeMultiPlasticityStress::singleStep(), TEST(), TEST_F(), ADComputeMultipleInelasticStress::updateQpStateSingleModel(), ComputeMultipleInelasticStressBase::updateQpStateSingleModel(), and MultiPlasticityRawComponentAssembler::yieldFunction().