Detailed Description

Definition at line 32 of file PhreeqcKernel.h.

#include <PhreeqcKernel.h>

Inheritance diagram for ChemistryLib::PhreeqcKernelData::PhreeqcKernel:

Collaboration diagram for ChemistryLib::PhreeqcKernelData::PhreeqcKernel:

Public Member Functions
	PhreeqcKernel (GlobalLinearSolver &linear_solver, std::size_t const num_chemical_systems, std::vector< std::pair< int, std::string >> const &process_id_to_component_name_map, std::string const &database, AqueousSolution aqueous_solution, std::unique_ptr< EquilibriumReactants > &&equilibrium_reactants, std::unique_ptr< Kinetics > &&kinetic_reactants, std::vector< ReactionRate > &&reaction_rates)

void	executeSpeciationCalculation (double const dt) override

void	setAqueousSolutions (std::vector< GlobalVector * > const &process_solutions)

void	callPhreeqc (std::vector< GlobalVector * > &process_solutions)

std::vector< GlobalVector * >	getIntPtProcessSolutions () const override

void	updateNodalProcessSolutions (std::vector< GlobalVector * > const &process_solutions, std::size_t const node_id)

Public Member Functions inherited from ChemistryLib::ChemicalSolverInterface
	ChemicalSolverInterface (GlobalLinearSolver &linear_solver_)

virtual void	initialize ()

virtual void	initializeChemicalSystemConcrete (std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)

virtual void	setChemicalSystemConcrete (std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const *, MaterialPropertyLib::VariableArray const &, ParameterLib::SpatialPosition const &, double const, double const)

virtual void	setAqueousSolutionsPrevFromDumpFile ()

virtual double	getConcentration (int const, GlobalIndexType const) const

virtual std::vector< std::string > const	getComponentList () const

virtual void	updateVolumeFractionPostReaction (GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const, double const, double const)

virtual void	updatePorosityPostReaction (GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)

virtual void	computeSecondaryVariable (std::size_t const, std::vector< GlobalIndexType > const &)

virtual	~ChemicalSolverInterface ()=default

Private Member Functions
void	initializePhreeqcGeneralSettings ()

void	tidyEquilibriumReactants (EquilibriumReactants const equilibrium_reactants)

void	loadDatabase (std::string const &database)

void	reinitializeRates ()

void	setConvergenceTolerance ()

void	configureOutputSettings ()

cxxISolution *	getOrCreateInitialAqueousSolution (cxxSolution &aqueous_solution)

bool	isHydrogen (char const *element) const

void	setTimeStepSize (double const dt)

void	reset (std::size_t const chemical_system_id)

Private Attributes
std::map< int, struct master * >	_process_id_to_master_map

std::unique_ptr< cxxISolution const >	_initial_aqueous_solution

std::unique_ptr< cxxSolution const >	_aqueous_solution

std::vector< ReactionRate > const	_reaction_rates

Additional Inherited Members
Public Attributes inherited from ChemistryLib::ChemicalSolverInterface
std::vector< GlobalIndexType >	chemical_system_index_map

GlobalLinearSolver &	linear_solver

Constructor & Destructor Documentation

◆ PhreeqcKernel()

ChemistryLib::PhreeqcKernelData::PhreeqcKernel::PhreeqcKernel	(	GlobalLinearSolver &	linear_solver,
		std::size_t const	num_chemical_systems,
		std::vector< std::pair< int, std::string >> const &	process_id_to_component_name_map,
		std::string const &	database,
		AqueousSolution	aqueous_solution,
		std::unique_ptr< EquilibriumReactants > &&	equilibrium_reactants,
		std::unique_ptr< Kinetics > &&	kinetic_reactants,
		std::vector< ReactionRate > &&	reaction_rates
	)

Definition at line 26 of file PhreeqcKernel.cpp.

     : ChemicalSolverInterface(linear_solver),
       _initial_aqueous_solution(aqueous_solution.getInitialAqueousSolution()),
       _aqueous_solution(aqueous_solution.castToBaseClassNoninitialized()),
       _reaction_rates(std::move(reaction_rates))
 {
     initializePhreeqcGeneralSettings();
  
     loadDatabase(database);
  
     // solution
     for (std::size_t chemical_system_id = 0;
          chemical_system_id < num_chemical_systems;
          ++chemical_system_id)
     {
         aqueous_solution.setChemicalSystemID(chemical_system_id);
         Rxn_solution_map.emplace(chemical_system_id,
                                  *aqueous_solution.castToBaseClass());
     }
     use.Set_solution_in(true);
  
     // equilibrium reactants
     if (equilibrium_reactants)
     {
         tidyEquilibriumReactants(*equilibrium_reactants);
  
         for (std::size_t chemical_system_id = 0;
              chemical_system_id < num_chemical_systems;
              ++chemical_system_id)
         {
             equilibrium_reactants->setChemicalSystemID(chemical_system_id);
             Rxn_pp_assemblage_map.emplace(
                 chemical_system_id, *equilibrium_reactants->castToBaseClass());
         }
         // explicitly release and delete equilibrium_reactants
         equilibrium_reactants.reset(nullptr);
  
         use.Set_pp_assemblage_in(true);
     }
  
     // kinetic reactants
     if (kinetic_reactants)
     {
         for (std::size_t chemical_system_id = 0;
              chemical_system_id < num_chemical_systems;
              ++chemical_system_id)
         {
             kinetic_reactants->setChemicalSystemID(chemical_system_id);
             Rxn_kinetics_map.emplace(chemical_system_id,
                                      *kinetic_reactants->castToBaseClass());
         }
         // explicitly release and delete kinetic_reactants
         kinetic_reactants.reset(nullptr);
  
         use.Set_kinetics_in(true);
     }
  
     // rates
     reinitializeRates();
  
     setConvergenceTolerance();
  
     configureOutputSettings();
  
     for (auto const& map_pair : process_id_to_component_name_map)
     {
         auto const transport_process_id = map_pair.first;
         auto const& transport_process_variable = map_pair.second;
  
         auto master_species =
             master_bsearch(transport_process_variable.c_str());
  
         _process_id_to_master_map[transport_process_id] = master_species;
     }
 }

References _process_id_to_master_map, ChemistryLib::PhreeqcKernelData::AqueousSolution::castToBaseClass(), configureOutputSettings(), initializePhreeqcGeneralSettings(), loadDatabase(), reinitializeRates(), ChemistryLib::PhreeqcKernelData::AqueousSolution::setChemicalSystemID(), setConvergenceTolerance(), and tidyEquilibriumReactants().

Member Function Documentation

◆ callPhreeqc()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::callPhreeqc ( std::vector< GlobalVector * > & process_solutions )

Definition at line 259 of file PhreeqcKernel.cpp.

 {
     std::size_t const num_chemical_systems = process_solutions[0]->size();
     for (std::size_t chemical_system_id = 0;
          chemical_system_id < num_chemical_systems;
          ++chemical_system_id)
     {
         Rxn_new_solution.insert(chemical_system_id);
         new_solution = 1;
         use.Set_n_solution_user(chemical_system_id);
  
         if (!Rxn_pp_assemblage_map.empty())
         {
             Rxn_new_pp_assemblage.insert(chemical_system_id);
             use.Set_pp_assemblage_in(true);
             use.Set_n_pp_assemblage_user(chemical_system_id);
         }
  
         if (!Rxn_kinetics_map.empty())
         {
             use.Set_kinetics_in(true);
             use.Set_n_kinetics_user(chemical_system_id);
         }
  
         initial_solutions(false);
  
         reactions();
  
         updateNodalProcessSolutions(process_solutions, chemical_system_id);
  
         reset(chemical_system_id);
     }
 }

References reset(), and updateNodalProcessSolutions().

Referenced by executeSpeciationCalculation().

◆ configureOutputSettings()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::configureOutputSettings ( )

inlineprivate

Definition at line 72 of file PhreeqcKernel.h.

72 { pr.all = false; }

Referenced by PhreeqcKernel().

◆ executeSpeciationCalculation()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::executeSpeciationCalculation ( double const dt )

overridevirtual

Implements ChemistryLib::ChemicalSolverInterface.

Definition at line 175 of file PhreeqcKernel.cpp.

 {
     std::vector<GlobalVector*> process_solutions;
  
     setAqueousSolutions(process_solutions);
  
     setTimeStepSize(dt);
  
     callPhreeqc(process_solutions);
 }

References callPhreeqc(), setAqueousSolutions(), and setTimeStepSize().

◆ getIntPtProcessSolutions()

std::vector<GlobalVector*> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::getIntPtProcessSolutions ( ) const

inlineoverridevirtual

Implements ChemistryLib::ChemicalSolverInterface.

Definition at line 51 of file PhreeqcKernel.h.

     {
         return {};
     }

◆ getOrCreateInitialAqueousSolution()

cxxISolution * ChemistryLib::PhreeqcKernelData::PhreeqcKernel::getOrCreateInitialAqueousSolution ( cxxSolution & aqueous_solution )

private

Definition at line 237 of file PhreeqcKernel.cpp.

 {
     if (!aqueous_solution.Get_initial_data())
     {
         aqueous_solution.Set_initial_data(_initial_aqueous_solution.get());
         aqueous_solution.Set_new_def(true);
     }
  
     return aqueous_solution.Get_initial_data();
 }

References _initial_aqueous_solution.

Referenced by setAqueousSolutions().

◆ initializePhreeqcGeneralSettings()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::initializePhreeqcGeneralSettings ( )

inlineprivate

Definition at line 61 of file PhreeqcKernel.h.

61 { do_initialize(); }

Referenced by PhreeqcKernel().

◆ isHydrogen()

bool ChemistryLib::PhreeqcKernelData::PhreeqcKernel::isHydrogen ( char const * element ) const

inlineprivate

Definition at line 77 of file PhreeqcKernel.h.

     {
         return strcmp(element, "H") == 0;
     }

Referenced by setAqueousSolutions(), and updateNodalProcessSolutions().

◆ loadDatabase()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::loadDatabase ( std::string const & database )

private

Definition at line 132 of file PhreeqcKernel.cpp.

 {
     std::ifstream in(database);
     if (!in)
     {
         OGS_FATAL("Unable to open database file '{:s}'.", database);
     }
     assert(phrq_io->get_istream() == nullptr);
     phrq_io->push_istream(&in, false);
     read_database();
 }

References OGS_FATAL.

Referenced by PhreeqcKernel().

◆ reinitializeRates()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::reinitializeRates ( )

private

Definition at line 144 of file PhreeqcKernel.cpp.

 {
     count_rates = _reaction_rates.size();
     rates = (struct rate*)realloc(
         rates, (std::size_t)(count_rates) * sizeof(struct rate));
     int rate_id = 0;
     for (auto const& reaction_rate : _reaction_rates)
     {
         rates[rate_id].name = reaction_rate.kinetic_reactant.data();
  
         // Commands' strings are freed by Phreeqc dtor.
         rates[rate_id].commands = static_cast<char*>(
             malloc(sizeof(char) * reaction_rate.commands().size() + 1));
         if (rates[rate_id].commands == nullptr)
         {
             OGS_FATAL("Could not allocate memory for rate[{:d}] commands.",
                       rate_id);
         }
         reaction_rate.commands().copy(rates[rate_id].commands,
                                       std::string::npos);
         rates[rate_id].commands[reaction_rate.commands().size()] = '\0';
  
         rates[rate_id].new_def = 1;
         rates[rate_id].linebase = nullptr;
         rates[rate_id].varbase = nullptr;
         rates[rate_id].loopbase = nullptr;
  
         ++rate_id;
     };
 }

References _reaction_rates, and OGS_FATAL.

Referenced by PhreeqcKernel().

◆ reset()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::reset ( std::size_t const chemical_system_id )

private

Definition at line 293 of file PhreeqcKernel.cpp.

 {
     // Clean up
     Rxn_new_solution.clear();
  
     // Solution
     {
         Rxn_solution_map[chemical_system_id] = *_aqueous_solution;
         Rxn_solution_map[chemical_system_id].Set_n_user_both(
             chemical_system_id);
         Rxn_solution_map[chemical_system_id].Set_pe(-s_eminus->la);
     }
  
     // Equilibrium reactants
     if (!Rxn_pp_assemblage_map.empty())
     {
         Rxn_new_pp_assemblage.clear();
         for (int j = 0; j < count_unknowns; j++)
         {
             if (x == nullptr || x[j]->type != PP)
                 continue;
  
             auto& phase_components = Rxn_pp_assemblage_map[chemical_system_id]
                                          .Get_pp_assemblage_comps();
             phase_components[x[j]->pp_assemblage_comp_name].Set_moles(
                 x[j]->moles);
         }
     }
  
     // Kinetics
     if (!Rxn_kinetics_map.empty())
     {
         Rxn_kinetics_map[chemical_system_id].Get_steps().clear();
     }
 }

References _aqueous_solution.

Referenced by callPhreeqc().

◆ setAqueousSolutions()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setAqueousSolutions ( std::vector< GlobalVector * > const & process_solutions )

Definition at line 186 of file PhreeqcKernel.cpp.

 {
     // Loop over chemical systems
     std::size_t const num_chemical_systems = process_solutions[0]->size();
     for (std::size_t chemical_system_id = 0;
          chemical_system_id < num_chemical_systems;
          ++chemical_system_id)
     {
         auto& aqueous_solution = Rxn_solution_map[chemical_system_id];
  
         auto initial_aqueous_solution =
             getOrCreateInitialAqueousSolution(aqueous_solution);
  
         auto& components = initial_aqueous_solution->Get_comps();
         // Loop over transport process id map to retrieve component
         // concentrations from process solutions
         for (auto const& map_pair : _process_id_to_master_map)
         {
             auto const transport_process_id = map_pair.first;
             auto const& master_species = map_pair.second;
  
             auto& transport_process_solution =
                 process_solutions[transport_process_id];
  
             auto& element_name = master_species->elt->name;
             auto const concentration =
                 transport_process_solution->get(chemical_system_id);
             if (isHydrogen(element_name))
             {
                 // Set pH value by hydrogen concentration.
                 double const pH = -std::log10(concentration);
                 aqueous_solution.Set_ph(pH);
  
                 {
                     auto hydrogen = components.find("H(1)");
                     if (hydrogen != components.end())
                     {
                         hydrogen->second.Set_input_conc(pH);
                     }
                 }
             }
             else
             {
                 // Set component concentrations.
                 components[element_name].Set_input_conc(concentration);
             }
         }
     }
 }

References _process_id_to_master_map, MaterialPropertyLib::concentration, getOrCreateInitialAqueousSolution(), isHydrogen(), and ChemistryLib::pH.

Referenced by executeSpeciationCalculation().

◆ setConvergenceTolerance()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setConvergenceTolerance ( )

inlineprivate

Definition at line 70 of file PhreeqcKernel.h.

70 { convergence_tolerance = 1e-12; }

Referenced by PhreeqcKernel().

◆ setTimeStepSize()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setTimeStepSize ( double const dt )

private

Definition at line 249 of file PhreeqcKernel.cpp.

 {
     // Loop over rxn kinetics map
     for (auto& map_pair : Rxn_kinetics_map)
     {
         auto& kinetics = map_pair.second;
         kinetics.Get_steps().push_back(dt);
     }
 }

Referenced by executeSpeciationCalculation().

◆ tidyEquilibriumReactants()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::tidyEquilibriumReactants ( EquilibriumReactants const equilibrium_reactants )

private

Definition at line 111 of file PhreeqcKernel.cpp.

 {
     // extract a part of function body from int
     // Phreeqc::tidy_pp_assemblage(void)
     count_elts = 0;
     double coef = 1.0;
     for (auto const& phase_component :
          equilibrium_reactants.getPhaseComponents())
     {
         int phase_id;
         struct phase* phase_component_ptr =
             phase_bsearch(phase_component.first.c_str(), &phase_id, FALSE);
         add_elt_list(phase_component_ptr->next_elt, coef);
     }
  
     cxxNameDouble nd = elt_list_NameDouble();
     const_cast<cxxPPassemblage*>(equilibrium_reactants.castToBaseClass())
         ->Set_eltList(nd);
 }

References ChemistryLib::PhreeqcKernelData::EquilibriumReactants::castToBaseClass(), and ChemistryLib::PhreeqcKernelData::EquilibriumReactants::getPhaseComponents().

Referenced by PhreeqcKernel().

◆ updateNodalProcessSolutions()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::updateNodalProcessSolutions	(	std::vector< GlobalVector * > const &	process_solutions,
		std::size_t const	node_id
	)

Definition at line 329 of file PhreeqcKernel.cpp.

 {
     for (auto const& map_pair : _process_id_to_master_map)
     {
         auto const transport_process_id = map_pair.first;
         auto const& master_species = map_pair.second;
  
         auto& transport_process_solution =
             process_solutions[transport_process_id];
  
         auto const& element_name = master_species->elt->name;
         if (isHydrogen(element_name))
         {
             // Update hydrogen concentration by pH value.
             auto const hydrogen_concentration = std::pow(10, s_hplus->la);
             transport_process_solution->set(node_id, hydrogen_concentration);
         }
         else
         {
             // Update solutions of component transport processes.
             auto const concentration =
                 master_species->primary
                     ? master_species->total_primary / mass_water_aq_x
                     : master_species->total / mass_water_aq_x;
             transport_process_solution->set(node_id, concentration);
         }
     }
 }

References _process_id_to_master_map, MaterialPropertyLib::concentration, and isHydrogen().

Referenced by callPhreeqc().

Member Data Documentation

◆ _aqueous_solution

std::unique_ptr<cxxSolution const> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_aqueous_solution

private

Definition at line 88 of file PhreeqcKernel.h.

Referenced by reset().

◆ _initial_aqueous_solution

std::unique_ptr<cxxISolution const> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_initial_aqueous_solution

private

Definition at line 87 of file PhreeqcKernel.h.

Referenced by getOrCreateInitialAqueousSolution().

◆ _process_id_to_master_map

std::map<int, struct master*> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_process_id_to_master_map

private

Definition at line 86 of file PhreeqcKernel.h.

Referenced by PhreeqcKernel(), setAqueousSolutions(), and updateNodalProcessSolutions().

◆ _reaction_rates

std::vector<ReactionRate> const ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_reaction_rates

private

Definition at line 89 of file PhreeqcKernel.h.

Referenced by reinitializeRates().

The documentation for this class was generated from the following files:

ChemistryLib/PhreeqcKernel.h
ChemistryLib/PhreeqcKernel.cpp

Detailed Description

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ PhreeqcKernel()

Member Function Documentation

◆ callPhreeqc()

◆ configureOutputSettings()

◆ executeSpeciationCalculation()

◆ getIntPtProcessSolutions()

◆ getOrCreateInitialAqueousSolution()

◆ initializePhreeqcGeneralSettings()

◆ isHydrogen()

◆ loadDatabase()

◆ reinitializeRates()

◆ reset()

◆ setAqueousSolutions()

◆ setConvergenceTolerance()

◆ setTimeStepSize()

◆ tidyEquilibriumReactants()

◆ updateNodalProcessSolutions()

Member Data Documentation

◆ _aqueous_solution

◆ _initial_aqueous_solution

◆ _process_id_to_master_map

◆ _reaction_rates