PhreeqcIO.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause

#pragma once
 
#include <memory>
 
#include "ChemicalSolverInterface.h"
#include "PhreeqcIOData/ClampingStats.h"
#include "PhreeqcIOData/Knobs.h"
#include "PhreeqcIOData/PhreeqcInstancePool.h"
 
namespace MeshLib
{
class Mesh;
}
 
namespace ChemistryLib
{
namespace PhreeqcIOData
{
struct ChemicalSystem;
struct ReactionRate;
struct Output;
struct Dump;
struct UserPunch;

class PhreeqcIO final : public ChemicalSolverInterface
{
public:
    PhreeqcIO(MeshLib::Mesh const& mesh,
              GlobalLinearSolver& linear_solver,
              std::string const& project_file_name,
              std::string&& database,
              std::unique_ptr<ChemicalSystem>&& chemical_system,
              std::vector<ReactionRate>&& reaction_rates,
              std::unique_ptr<UserPunch>&& user_punch,
              std::unique_ptr<Output>&& output,
              std::unique_ptr<Dump>&& dump,
              Knobs&& knobs,
              bool use_stream_mode,
              int num_chemistry_threads,
              double concentration_warning_threshold);
 
    ~PhreeqcIO();
 
    void initialize() override;
 
    void initializeChemicalSystemConcrete(
        std::vector<double> const& concentrations,
        GlobalIndexType const& chemical_system_id,
        MaterialPropertyLib::Medium const& medium,
        ParameterLib::SpatialPosition const& pos,
        double const t) override;
 
    void setChemicalSystemConcrete(
        std::vector<double> const& concentrations,
        GlobalIndexType const& chemical_system_id,
        MaterialPropertyLib::Medium const* medium,
        MaterialPropertyLib::VariableArray const& vars,
        ParameterLib::SpatialPosition const& pos, double const t,
        double const dt) override;
 
    void setAqueousSolutionsPrevFromDumpFile() override;
 
    void executeSpeciationCalculation(double const dt) override;
 
    double getConcentration(
        int const component_id,
        GlobalIndexType const chemical_system_id) const override;
 
    friend std::ostream& operator<<(std::ostream& os,
                                    PhreeqcIO const& phreeqc_io);
 
    friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
 
    void updateVolumeFractionPostReaction(
        GlobalIndexType const& chemical_system_id,
        MaterialPropertyLib::Medium const& medium,
        ParameterLib::SpatialPosition const& pos, double const porosity,
        double const t, double const dt) override;
 
    void updatePorosityPostReaction(GlobalIndexType const& chemical_system_id,
                                    MaterialPropertyLib::Medium const& medium,
                                    double& porosity) override;
 
    void computeSecondaryVariable(
        std::size_t const ele_id,
        std::vector<GlobalIndexType> const& chemical_system_indices) override;
 
    std::vector<std::string> const getComponentList() const override;
 
    std::string const _phreeqc_input_file;
 
private:
    void writeInputsToFile(double const dt);

    void writeInputHeader(std::ostream& os) const;

    void writeSystemBlock(std::ostream& os,
                          std::size_t chemical_system_id,
                          std::size_t solution_id,
                          std::size_t prev_solution_id,
                          double dt) const;

    void updateSystemFromOutputLine(std::string_view line,
                                    std::size_t chemical_system_id);
 
    void setAqueousSolutionsPrevFromDumpString(std::string_view dump_content);
 
    void callPhreeqc() const;
 
    void readOutputsFromFile();
 
    void executeSpeciationCalculationParallel(double const dt);
 
    std::string generateInputForSystem(std::size_t chemical_system_id,
                                       double const dt) const;
 
    void parseOutputForSystem(std::string_view output_content,
                              std::size_t chemical_system_id);
 
    void updateChemicalSystemFromOutput(
        std::vector<double> const& accepted_items,
        std::size_t chemical_system_id);
 
    PhreeqcIO& operator<<(double const dt)
    {
        _dt = dt;
        return *this;
    }
 
    // Member variables are ordered by size (largest first) to minimize padding.

    ClampingStats _clamping_totals;
    std::string const _database;
    Knobs const _knobs;
    std::vector<ReactionRate> const _reaction_rates;
    std::unique_ptr<ChemicalSystem> _chemical_system;
    std::unique_ptr<UserPunch> _user_punch;
    std::unique_ptr<Output> const _output;
    std::unique_ptr<Dump> const _dump;
    std::unique_ptr<PhreeqcInstancePool> instance_pool_;
    double _dt = std::numeric_limits<double>::quiet_NaN();
    std::size_t _num_chemical_systems = 0;
    double _concentration_warning_threshold;
    int phreeqc_instance_id = -1;
    int num_chemistry_threads_;
    bool _use_stream_mode;
};
}  // namespace PhreeqcIOData
}  // namespace ChemistryLib