PertinentGeochemicalSystem.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "GeochemicalDatabaseReader.h"
#include "GeochemistryKineticRateCalculator.h"
#include <unordered_map>
#include "DenseMatrix.h"
#include <libmesh/dense_vector.h>
 
struct SurfaceComplexationInfo
{
  SurfaceComplexationInfo() {}
 
  bool operator==(const SurfaceComplexationInfo & rhs) const
  {
    return surface_area == rhs.surface_area && sorption_sites == rhs.sorption_sites;
  };
 
  Real surface_area;
  std::map<std::string, Real> sorption_sites;
};
 
struct KineticRateDefinition
{
  KineticRateDefinition(unsigned kinetic_species_index,
                        const std::vector<Real> & promoting_indices,
                        const std::vector<Real> & promoting_monod_indices,
                        const std::vector<Real> & promoting_half_saturation,
                        unsigned progeny_index,
                        const KineticRateUserDescription & description)
    : kinetic_species_index(kinetic_species_index),
      promoting_indices(promoting_indices),
      promoting_monod_indices(promoting_monod_indices),
      promoting_half_saturation(promoting_half_saturation),
      progeny_index(progeny_index),
      description(description)
  {
  }
 
  bool operator==(const KineticRateDefinition & rhs) const
  {
    return (kinetic_species_index == rhs.kinetic_species_index) &&
           (promoting_indices == rhs.promoting_indices) &&
           (promoting_monod_indices == rhs.promoting_monod_indices) &&
           (promoting_half_saturation == rhs.promoting_half_saturation) &&
           (progeny_index == rhs.progeny_index) && (description == rhs.description);
  };
 
  unsigned kinetic_species_index;
  std::vector<Real> promoting_indices;
  std::vector<Real> promoting_monod_indices;
  std::vector<Real> promoting_half_saturation;
  unsigned progeny_index;
  KineticRateUserDescription description;
};
 
struct ModelGeochemicalDatabase
{
  ModelGeochemicalDatabase(const GeochemicalDatabaseReader & db)
    : original_database(&db), swap_to_original_basis(DenseMatrix<Real>())
  {
  }
 
  bool operator==(const ModelGeochemicalDatabase & rhs) const
  {
    return (original_database == rhs.original_database) &&
           (basis_species_index == rhs.basis_species_index) &&
           (basis_species_name == rhs.basis_species_name) &&
           (basis_species_mineral == rhs.basis_species_mineral) &&
           (basis_species_gas == rhs.basis_species_gas) &&
           (basis_species_transported == rhs.basis_species_transported) &&
           (basis_species_charge == rhs.basis_species_charge) &&
           (basis_species_radius == rhs.basis_species_radius) &&
           (basis_species_molecular_weight == rhs.basis_species_molecular_weight) &&
           (basis_species_molecular_volume == rhs.basis_species_molecular_volume) &&
           (eqm_species_index == rhs.eqm_species_index) &&
           (eqm_species_name == rhs.eqm_species_name) &&
           (eqm_species_mineral == rhs.eqm_species_mineral) &&
           (eqm_species_gas == rhs.eqm_species_gas) &&
           (eqm_species_transported == rhs.eqm_species_transported) &&
           (eqm_species_charge == rhs.eqm_species_charge) &&
           (eqm_species_radius == rhs.eqm_species_radius) &&
           (eqm_species_molecular_weight == rhs.eqm_species_molecular_weight) &&
           (eqm_species_molecular_volume == rhs.eqm_species_molecular_volume) &&
           (eqm_stoichiometry == rhs.eqm_stoichiometry) && (eqm_log10K == rhs.eqm_log10K) &&
           (surface_sorption_name == rhs.surface_sorption_name) &&
           (surface_sorption_area == rhs.surface_sorption_area) &&
           (surface_sorption_related == rhs.surface_sorption_related) &&
           (surface_sorption_number == rhs.surface_sorption_number) &&
           (redox_lhs == rhs.redox_lhs) && (redox_stoichiometry == rhs.redox_stoichiometry) &&
           (redox_log10K == rhs.redox_log10K) &&
           (surface_complexation_info == rhs.surface_complexation_info) &&
           (gas_chi == rhs.gas_chi) && (kin_species_index == rhs.kin_species_index) &&
           (kin_species_name == rhs.kin_species_name) &&
           (kin_species_mineral == rhs.kin_species_mineral) &&
           (kin_species_transported == rhs.kin_species_transported) &&
           (kin_species_charge == rhs.kin_species_charge) &&
           (kin_species_molecular_weight == rhs.kin_species_molecular_weight) &&
           (kin_species_molecular_volume == rhs.kin_species_molecular_volume) &&
           (kin_log10K == rhs.kin_log10K) && (kin_stoichiometry == rhs.kin_stoichiometry) &&
           (kin_rate == rhs.kin_rate) &&
           (have_swapped_out_of_basis == rhs.have_swapped_out_of_basis) &&
           (have_swapped_into_basis == rhs.have_swapped_into_basis) &&
           (swap_to_original_basis == rhs.swap_to_original_basis);
  };
 
  const GeochemicalDatabaseReader * original_database;
 
  std::unordered_map<std::string, unsigned> basis_species_index;
 
  std::vector<std::string> basis_species_name;
 
  std::vector<bool> basis_species_mineral;
 
  std::vector<bool> basis_species_gas;
 
  std::vector<bool> basis_species_transported;
 
  std::vector<Real> basis_species_charge;
 
  std::vector<Real> basis_species_radius;
 
  std::vector<Real> basis_species_molecular_weight;
 
  std::vector<Real> basis_species_molecular_volume;
 
  std::unordered_map<std::string, unsigned> eqm_species_index;
 
  std::vector<std::string> eqm_species_name;
 
  std::vector<bool> eqm_species_mineral;
 
  std::vector<bool> eqm_species_gas;
 
  std::vector<bool> eqm_species_transported;
 
  std::vector<Real> eqm_species_charge;
 
  std::vector<Real> eqm_species_radius;
 
  std::vector<Real> eqm_species_molecular_weight;
 
  std::vector<Real> eqm_species_molecular_volume;
 
  DenseMatrix<Real> eqm_stoichiometry;
 
  DenseMatrix<Real> eqm_log10K;
 
  std::vector<std::string> surface_sorption_name;
 
  std::vector<Real> surface_sorption_area;
 
  std::vector<bool> surface_sorption_related;
 
  std::vector<unsigned> surface_sorption_number;
 
  std::string redox_lhs;
 
  DenseMatrix<Real> redox_stoichiometry;
 
  DenseMatrix<Real> redox_log10K;
 
  std::unordered_map<std::string, SurfaceComplexationInfo> surface_complexation_info;
 
  std::unordered_map<std::string, std::vector<Real>> gas_chi;
 
  std::unordered_map<std::string, unsigned> kin_species_index;
 
  std::vector<std::string> kin_species_name;
 
  std::vector<bool> kin_species_mineral;
 
  std::vector<bool> kin_species_transported;
 
  std::vector<Real> kin_species_charge;
 
  std::vector<Real> kin_species_molecular_weight;
 
  std::vector<Real> kin_species_molecular_volume;
 
  DenseMatrix<Real> kin_log10K;
 
  DenseMatrix<Real> kin_stoichiometry;
 
  std::vector<KineticRateDefinition> kin_rate;
 
  std::vector<unsigned> have_swapped_out_of_basis;
 
  std::vector<unsigned> have_swapped_into_basis;
 
  DenseMatrix<Real> swap_to_original_basis;
};
 
class PertinentGeochemicalSystem
{
public:
  PertinentGeochemicalSystem(const GeochemicalDatabaseReader & db,
                             const std::vector<std::string> & basis_species,
                             const std::vector<std::string> & minerals,
                             const std::vector<std::string> & gases,
                             const std::vector<std::string> & kinetic_minerals,
                             const std::vector<std::string> & kinetic_redox,
                             const std::vector<std::string> & kinetic_surface_species,
                             const std::string & redox_ox,
                             const std::string & redox_e);
 
  const ModelGeochemicalDatabase & modelGeochemicalDatabase() const;
 
  void addKineticRate(const KineticRateUserDescription & description);
 
  unsigned getIndexOfOriginalBasisSpecies(const std::string & name) const;
 
  std::vector<std::string> originalBasisNames() const;
 
private:
  GeochemicalDatabaseReader _db;
 
  std::unordered_map<std::string, unsigned> _basis_index;
 
  std::vector<GeochemistryBasisSpecies> _basis_info;
 
  std::unordered_map<std::string, unsigned> _mineral_index;
 
  std::vector<GeochemistryMineralSpecies> _mineral_info;
 
  std::unordered_map<std::string, unsigned> _gas_index;
 
  std::vector<GeochemistryGasSpecies> _gas_info;
 
  std::unordered_map<std::string, unsigned> _kinetic_mineral_index;
 
  std::vector<GeochemistryMineralSpecies> _kinetic_mineral_info;
 
  std::unordered_map<std::string, unsigned> _kinetic_redox_index;
 
  std::vector<GeochemistryRedoxSpecies> _kinetic_redox_info;
 
  std::unordered_map<std::string, unsigned> _kinetic_surface_index;
 
  std::vector<GeochemistrySurfaceSpecies> _kinetic_surface_info;
 
  std::unordered_map<std::string, unsigned> _secondary_index;
 
  std::vector<GeochemistryEquilibriumSpecies> _secondary_info;
 
  const std::string _redox_ox;
 
  const std::string _redox_e;
 
  ModelGeochemicalDatabase _model;
 
  void buildBasis(const std::vector<std::string> & basis_species);
 
  void buildMinerals(const std::vector<std::string> & minerals);
 
  void buildAllMinerals(const std::vector<std::string> & minerals);
 
  void buildGases(const std::vector<std::string> & gases);
 
  void buildKineticMinerals(const std::vector<std::string> & kinetic_minerals);
 
  void buildKineticRedox(const std::vector<std::string> & kinetic_redox);
 
  void buildKineticSurface(const std::vector<std::string> & kinetic_surface);
 
  void buildSecondarySpecies();
 
  bool checkRedoxe();
 
  void checkMinerals(const std::vector<GeochemistryMineralSpecies> & mineral_info) const;
 
  void checkGases() const;
 
  void checkKineticRedox() const;
 
  void checkKineticSurfaceSpecies() const;
 
  void createModel();
 
  void buildRedoxeInfo(std::vector<Real> & redox_e_stoichiometry,
                       std::vector<Real> & redox_e_log10K);
};