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);
};