diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index 7212fda32ab4c663ea10b9492a85469358048a1a..11f677e6121ea016ae1645a21061241a5b3681f6 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -17,17 +17,13 @@
#include "Common/CreateReactionRate.h"
#include "MeshLib/Mesh.h"
#include "PhreeqcIO.h"
-#include "PhreeqcIOData/AqueousSolution.h"
-#include "PhreeqcIOData/CreateAqueousSolution.h"
-#include "PhreeqcIOData/CreateEquilibriumReactants.h"
-#include "PhreeqcIOData/CreateKineticReactant.h"
+#include "PhreeqcIOData/ChemicalSystem.h"
+#include "PhreeqcIOData/CreateChemicalSystem.h"
#include "PhreeqcIOData/CreateKnobs.h"
#include "PhreeqcIOData/CreateOutput.h"
#include "PhreeqcIOData/CreateSurface.h"
#include "PhreeqcIOData/CreateUserPunch.h"
#include "PhreeqcIOData/Dump.h"
-#include "PhreeqcIOData/EquilibriumReactant.h"
-#include "PhreeqcIOData/KineticReactant.h"
#include "PhreeqcIOData/Knobs.h"
#include "PhreeqcIOData/ReactionRate.h"
#include "PhreeqcIOData/Surface.h"
@@ -87,32 +83,18 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
auto path_to_database = parseDatabasePath(config);
- // solution
- auto aqueous_solution = PhreeqcIOData::createAqueousSolution(
- //! \ogs_file_param{prj__chemical_system__solution}
- config.getConfigSubtree("solution"));
-
- // kinetic reactants
+ // chemical system
auto chemical_system_map =
*mesh.getProperties().template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);
-
- auto kinetic_reactants = PhreeqcIOData::createKineticReactants(
- //! \ogs_file_param{prj__chemical_system__kinetic_reactants}
- config.getConfigSubtreeOptional("kinetic_reactants"), *meshes[0],
- chemical_system_map);
+ auto chemical_system = PhreeqcIOData::createChemicalSystem(
+ config, *meshes[0], chemical_system_map);
// rates
auto reaction_rates = createReactionRates<PhreeqcIOData::ReactionRate>(
//! \ogs_file_param{prj__chemical_system__rates}
config.getConfigSubtreeOptional("rates"));
- // equilibrium reactants
- auto equilibrium_reactants = PhreeqcIOData::createEquilibriumReactants(
- //! \ogs_file_param{prj__chemical_system__equilibrium_reactants}
- config.getConfigSubtreeOptional("equilibrium_reactants"), *meshes[0],
- chemical_system_map);
-
// surface
auto surface = PhreeqcIOData::createSurface(
//! \ogs_file_param{prj__chemical_system__surface}
@@ -138,22 +120,15 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
config.getConfigSubtreeOptional("user_punch"), *meshes[0]);
// output
- auto const& components = aqueous_solution.components;
auto const use_high_precision =
//! \ogs_file_param{prj__chemical_system__use_high_precision}
config.getConfigParameter<bool>("use_high_precision", true);
auto output = PhreeqcIOData::createOutput(
- components, equilibrium_reactants, kinetic_reactants, user_punch,
- use_high_precision, project_file_name);
-
- auto const num_chemical_systems = mesh.getNumberOfBaseNodes();
- std::vector<PhreeqcIOData::AqueousSolution> aqueous_solutions(
- num_chemical_systems, aqueous_solution);
+ *chemical_system, user_punch, use_high_precision, project_file_name);
return std::make_unique<PhreeqcIOData::PhreeqcIO>(
std::move(project_file_name), *meshes[mesh.getID()],
- std::move(path_to_database), std::move(aqueous_solutions),
- std::move(equilibrium_reactants), std::move(kinetic_reactants),
+ std::move(path_to_database), std::move(chemical_system),
std::move(reaction_rates), std::move(surface), std::move(user_punch),
std::move(output), std::move(dump), std::move(knobs));
}
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index cb6088a50e719a9712214df1d00c986d5412ac16..54ec132c52c9a7c9708803fae7fe39ee9d75fea9 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -18,8 +18,10 @@
#include "BaseLib/Algorithm.h"
#include "BaseLib/ConfigTreeUtil.h"
+#include "MathLib/LinAlg/Eigen/EigenVector.h"
#include "MeshLib/Mesh.h"
#include "PhreeqcIOData/AqueousSolution.h"
+#include "PhreeqcIOData/ChemicalSystem.h"
#include "PhreeqcIOData/Dump.h"
#include "PhreeqcIOData/EquilibriumReactant.h"
#include "PhreeqcIOData/KineticReactant.h"
@@ -48,9 +50,7 @@ std::ostream& operator<<(std::ostream& os,
PhreeqcIO::PhreeqcIO(std::string const project_file_name,
MeshLib::Mesh const& mesh,
std::string&& database,
- std::vector<AqueousSolution>&& aqueous_solutions,
- std::vector<EquilibriumReactant>&& equilibrium_reactants,
- std::vector<KineticReactant>&& kinetic_reactants,
+ std::unique_ptr<ChemicalSystem>&& chemical_system,
std::vector<ReactionRate>&& reaction_rates,
std::vector<SurfaceSite>&& surface,
std::unique_ptr<UserPunch>&& user_punch,
@@ -60,9 +60,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name,
: _phreeqc_input_file(project_file_name + "_phreeqc.inp"),
_mesh(mesh),
_database(std::move(database)),
- _aqueous_solutions(std::move(aqueous_solutions)),
- _equilibrium_reactants(std::move(equilibrium_reactants)),
- _kinetic_reactants(std::move(kinetic_reactants)),
+ _chemical_system(std::move(chemical_system)),
_reaction_rates(std::move(reaction_rates)),
_surface(std::move(surface)),
_user_punch(std::move(user_punch)),
@@ -147,16 +145,14 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);
// Loop over chemical systems
+ auto& aqueous_solution = _chemical_system->aqueous_solution;
+ auto& components = aqueous_solution->components;
for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
{
auto const global_id = chemical_system_map[local_id];
- // Get chemical compostion of solution in a particular chemical system
- auto& aqueous_solution = _aqueous_solutions[local_id];
- auto& components = aqueous_solution.components;
// Loop over transport process id map to retrieve component
// concentrations from process solutions or to update process solutions
// after chemical calculation by Phreeqc
-
for (unsigned component_id = 0; component_id < components.size();
++component_id)
{
@@ -167,13 +163,13 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
{
case Status::SettingAqueousSolutions:
// Set component concentrations.
- component.amount =
- transport_process_solution->get(global_id);
+ component.amount->set(
+ local_id, (*transport_process_solution)[global_id]);
break;
case Status::UpdatingProcessSolutions:
// Update solutions of component transport processes.
- transport_process_solution->set(global_id,
- component.amount);
+ transport_process_solution->set(
+ global_id, (*component.amount)[local_id]);
break;
}
}
@@ -183,15 +179,16 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
case Status::SettingAqueousSolutions:
{
// Set pH value by hydrogen concentration.
- aqueous_solution.pH =
- -std::log10(process_solutions.back()->get(global_id));
+ aqueous_solution->pH->set(
+ local_id,
+ -std::log10((*process_solutions.back())[global_id]));
break;
}
case Status::UpdatingProcessSolutions:
{
// Update hydrogen concentration by pH value.
auto hydrogen_concentration =
- std::pow(10, -aqueous_solution.pH);
+ std::pow(10, -(*aqueous_solution->pH)[local_id]);
process_solutions.back()->set(global_id,
hydrogen_concentration);
break;
@@ -277,10 +274,8 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
{
auto const global_id = chemical_system_map[local_id];
- auto const& aqueous_solution =
- phreeqc_io._aqueous_solutions[local_id];
os << "SOLUTION " << global_id + 1 << "\n";
- os << aqueous_solution << "\n";
+ phreeqc_io._chemical_system->aqueous_solution->print(os, local_id);
auto const& dump = phreeqc_io._dump;
if (dump)
@@ -297,7 +292,8 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
os << "USE solution " << global_id + 1 << "\n\n";
- auto const& equilibrium_reactants = phreeqc_io._equilibrium_reactants;
+ auto const& equilibrium_reactants =
+ phreeqc_io._chemical_system->equilibrium_reactants;
if (!equilibrium_reactants.empty())
{
os << "EQUILIBRIUM_PHASES " << global_id + 1 << "\n";
@@ -308,7 +304,8 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
os << "\n";
}
- auto const& kinetic_reactants = phreeqc_io._kinetic_reactants;
+ auto const& kinetic_reactants =
+ phreeqc_io._chemical_system->kinetic_reactants;
if (!kinetic_reactants.empty())
{
os << "KINETICS " << global_id + 1 << "\n";
@@ -461,11 +458,12 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
}
assert(accepted_items.size() == output.accepted_items.size());
- auto& aqueous_solution =
- phreeqc_io._aqueous_solutions[local_id];
- auto& components = aqueous_solution.components;
- auto& equilibrium_reactants = phreeqc_io._equilibrium_reactants;
- auto& kinetic_reactants = phreeqc_io._kinetic_reactants;
+ auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
+ auto& components = aqueous_solution->components;
+ auto& equilibrium_reactants =
+ phreeqc_io._chemical_system->equilibrium_reactants;
+ auto& kinetic_reactants =
+ phreeqc_io._chemical_system->kinetic_reactants;
auto& user_punch = phreeqc_io._user_punch;
for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
++item_id)
@@ -482,13 +480,15 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
case ItemType::pH:
{
// Update pH value
- aqueous_solution.pH = accepted_items[item_id];
+ aqueous_solution->pH->set(local_id,
+ accepted_items[item_id]);
break;
}
case ItemType::pe:
{
// Update pe value
- aqueous_solution.pe = accepted_items[item_id];
+ (*aqueous_solution->pe)[global_id] =
+ accepted_items[item_id];
break;
}
case ItemType::Component:
@@ -497,7 +497,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
auto& component = BaseLib::findElementOrError(
components.begin(), components.end(), compare_by_name,
"Could not find component '" + item_name + "'.");
- component.amount = accepted_items[item_id];
+ component.amount->set(local_id, accepted_items[item_id]);
break;
}
case ItemType::EquilibriumReactant:
@@ -547,7 +547,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
std::vector<std::string> const PhreeqcIO::getComponentList() const
{
std::vector<std::string> component_names;
- auto const& components = _aqueous_solutions.front().components;
+ auto const& components = _chemical_system->aqueous_solution->components;
std::transform(components.begin(), components.end(),
std::back_inserter(component_names),
[](auto const& c) { return c.name; });
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 1835c743aa7eb2cbf69c9e35801f936ee9a35f0c..4c59ca563b295b92f38f76e7ebc2aa1066177baf 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -24,9 +24,7 @@ namespace ChemistryLib
{
namespace PhreeqcIOData
{
-struct AqueousSolution;
-struct EquilibriumReactant;
-struct KineticReactant;
+struct ChemicalSystem;
struct ReactionRate;
struct Output;
struct SurfaceSite;
@@ -45,9 +43,7 @@ public:
PhreeqcIO(std::string const project_file_name,
MeshLib::Mesh const& mesh,
std::string&& database,
- std::vector<AqueousSolution>&& aqueous_solutions,
- std::vector<EquilibriumReactant>&& equilibrium_reactants,
- std::vector<KineticReactant>&& kinetic_reactants,
+ std::unique_ptr<ChemicalSystem>&& chemical_system,
std::vector<ReactionRate>&& reaction_rates,
std::vector<SurfaceSite>&& surface,
std::unique_ptr<UserPunch>&& user_punch,
@@ -92,9 +88,7 @@ private:
MeshLib::Mesh const& _mesh;
std::string const _database;
- std::vector<AqueousSolution> _aqueous_solutions;
- std::vector<EquilibriumReactant> _equilibrium_reactants;
- std::vector<KineticReactant> _kinetic_reactants;
+ std::unique_ptr<ChemicalSystem> _chemical_system;
std::vector<ReactionRate> const _reaction_rates;
std::vector<SurfaceSite> const _surface;
std::unique_ptr<UserPunch> _user_punch;
diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
index fe0a5d0f18c6efa705f14093c12e7effba782c08..5d80ad5023415b69fe3935d07b58db8ad5e2f25f 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
@@ -12,44 +12,46 @@
#include "AqueousSolution.h"
#include "ChemistryLib/Common/ChargeBalance.h"
+#include "MeshLib/PropertyVector.h"
namespace ChemistryLib
{
namespace PhreeqcIOData
{
-std::ostream& operator<<(std::ostream& os,
- AqueousSolution const& aqueous_solution)
+void AqueousSolution::print(std::ostream& os,
+ std::size_t const chemical_system_id) const
{
- os << "temp " << aqueous_solution.temperature << "\n";
+ os << "temp " << temperature << "\n";
- os << "pressure " << aqueous_solution.pressure << "\n";
+ os << "pressure " << pressure << "\n";
- switch (aqueous_solution.charge_balance)
+ switch (charge_balance)
{
case ChargeBalance::pH:
- os << "pH " << aqueous_solution.pH << " charge"
+ os << "pH " << (*pH)[chemical_system_id] << " charge"
<< "\n";
- os << "pe " << aqueous_solution.pe << "\n";
+ os << "pe " << (*pe)[chemical_system_id] << "\n";
break;
case ChargeBalance::pe:
- os << "pH " << aqueous_solution.pH << "\n";
- os << "pe " << aqueous_solution.pe << " charge"
+ os << "pH " << (*pH)[chemical_system_id] << "\n";
+ os << "pe " << (*pe)[chemical_system_id] << " charge"
<< "\n";
break;
case ChargeBalance::Unspecified:
- os << "pH " << aqueous_solution.pH << "\n";
- os << "pe " << aqueous_solution.pe << "\n";
+ os << "pH " << (*pH)[chemical_system_id] << "\n";
+ os << "pe " << (*pe)[chemical_system_id] << "\n";
break;
}
os << "units mol/kgw\n";
- for (auto const& component : aqueous_solution.components)
+ for (auto const& component : components)
{
- os << component.name << " " << component.amount << "\n";
+ os << component.name << " " << (*component.amount)[chemical_system_id]
+ << "\n";
}
- return os;
+ os << "\n\n";
}
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.h b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
index 715a58a3d62bf1c68e66aea2244b76afbe56758d..e7c8d7f14f720e2030d0c934b8b86ae868ced9b6 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.h
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
@@ -11,11 +11,20 @@
#pragma once
#include <iosfwd>
+#include <memory>
#include <string>
#include <vector>
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "Output.h"
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
namespace ChemistryLib
{
enum class ChargeBalance;
@@ -24,35 +33,42 @@ namespace PhreeqcIOData
{
struct Component
{
- explicit Component(std::string name_) : name(std::move(name_)) {}
+ explicit Component(std::string name_,
+ std::size_t const num_chemical_systems_)
+ : name(std::move(name_)),
+ amount(MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ num_chemical_systems_))
+ {
+ }
std::string const name;
- double amount = std::numeric_limits<double>::quiet_NaN();
+ std::unique_ptr<GlobalVector> amount;
static const ItemType item_type = ItemType::Component;
};
struct AqueousSolution
{
- AqueousSolution(double temperature_,
- double pressure_,
- double pe_,
+ AqueousSolution(double temperature_, double pressure_,
+ MeshLib::PropertyVector<double>* pe_,
std::vector<Component>&& components_,
- ChargeBalance charge_balance_)
+ ChargeBalance charge_balance_,
+ std::size_t const num_chemical_systems_)
: temperature(temperature_),
pressure(pressure_),
+ pH(MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ num_chemical_systems_)),
pe(pe_),
components(std::move(components_)),
charge_balance(charge_balance_)
{
}
- friend std::ostream& operator<<(std::ostream& os,
- AqueousSolution const& aqueous_solution);
+ void print(std::ostream& os, std::size_t const chemical_system_id) const;
- double temperature;
- double pressure;
- double pH = std::numeric_limits<double>::quiet_NaN();
- double pe;
+ double const temperature;
+ double const pressure;
+ std::unique_ptr<GlobalVector> pH;
+ MeshLib::PropertyVector<double>* pe;
std::vector<Component> components;
ChargeBalance const charge_balance;
};
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.h b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..596c4ad2f63f5e566da37b9df8647bfe0f0c982e
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
@@ -0,0 +1,45 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "AqueousSolution.h"
+#include "EquilibriumReactant.h"
+#include "KineticReactant.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+struct ChemicalSystem
+{
+ ChemicalSystem(std::unique_ptr<AqueousSolution>&& aqueous_solution_,
+ std::vector<KineticReactant>&& kinetic_reactants_,
+ std::vector<EquilibriumReactant>&& equilibrium_reactants_)
+ : aqueous_solution(std::move(aqueous_solution_)),
+ kinetic_reactants(std::move(kinetic_reactants_)),
+ equilibrium_reactants(std::move(equilibrium_reactants_))
+ {
+ }
+
+ std::unique_ptr<AqueousSolution> aqueous_solution;
+ std::vector<KineticReactant> kinetic_reactants;
+ std::vector<EquilibriumReactant> equilibrium_reactants;
+};
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
index 70a7e2c6157b7689adb68b372e72d35912731303..91ba7f9015768c1bec1fde5aea1ecebea13fbb72 100644
--- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
@@ -13,12 +13,16 @@
#include "BaseLib/ConfigTree.h"
#include "ChemistryLib/Common/CreateChargeBalance.h"
#include "CreateSolutionComponent.h"
+#include "MeshLib/Mesh.h"
namespace ChemistryLib
{
namespace PhreeqcIOData
{
-AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config)
+std::unique_ptr<AqueousSolution> createAqueousSolution(
+ BaseLib::ConfigTree const& config,
+ MeshLib::Mesh const& mesh,
+ MeshLib::PropertyVector<std::size_t> const& chemical_system_map)
{
//! \ogs_file_param{prj__chemical_system__solution__temperature}
auto const temperature = config.getConfigParameter<double>("temperature");
@@ -27,13 +31,31 @@ AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config)
auto const pressure = config.getConfigParameter<double>("pressure");
//! \ogs_file_param{prj__chemical_system__solution__pe}
- auto const pe = config.getConfigParameter<double>("pe");
+ auto const pe0 = config.getConfigParameter<double>("pe");
- auto components = createSolutionComponents(config);
+ auto pe = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "pe", MeshLib::MeshItemType::Node, 1);
+
+ std::fill(std::begin(*pe),
+ std::end(*pe),
+ std::numeric_limits<double>::quiet_NaN());
+
+ std::for_each(
+ chemical_system_map.begin(),
+ chemical_system_map.end(),
+ [&pe, pe0](auto const& global_id) { (*pe)[global_id] = pe0; });
+
+ auto components =
+ createSolutionComponents(config, mesh.getNumberOfBaseNodes());
auto charge_balance = createChargeBalance(config);
- return {temperature, pressure, pe, std::move(components), charge_balance};
+ return std::make_unique<AqueousSolution>(temperature,
+ pressure,
+ pe,
+ std::move(components),
+ charge_balance,
+ mesh.getNumberOfBaseNodes());
}
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
index 8e0d680f276eca3a44bf227a0a6c7c9576db7a21..d83f595b92433cc704bbb7390dc76250b60de46e 100644
--- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
+++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
@@ -10,17 +10,30 @@
#pragma once
+#include <memory>
+
namespace BaseLib
{
class ConfigTree;
}
+namespace MeshLib
+{
+class Mesh;
+
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+} // namespace MeshLib
+
namespace ChemistryLib
{
namespace PhreeqcIOData
{
struct AqueousSolution;
-AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config);
+std::unique_ptr<AqueousSolution> createAqueousSolution(
+ BaseLib::ConfigTree const& config,
+ MeshLib::Mesh const& mesh,
+ MeshLib::PropertyVector<std::size_t> const& chemical_system_map);
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8cc8ee7612f6b1ee6f821542f936ef8f8c3ffa8
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.cpp
@@ -0,0 +1,50 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateChemicalSystem.h"
+#include "BaseLib/ConfigTree.h"
+#include "ChemicalSystem.h"
+#include "CreateAqueousSolution.h"
+#include "CreateEquilibriumReactants.h"
+#include "CreateKineticReactant.h"
+#include "EquilibriumReactant.h"
+#include "KineticReactant.h"
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+std::unique_ptr<ChemicalSystem> createChemicalSystem(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh,
+ MeshLib::PropertyVector<std::size_t> const& chemical_system_map)
+{
+ // solution
+ auto aqueous_solution = createAqueousSolution(
+ //! \ogs_file_param{prj__chemical_system__solution}
+ config.getConfigSubtree("solution"), mesh, chemical_system_map);
+
+ // kinetic reactants
+ auto kinetic_reactants = createKineticReactants(
+ //! \ogs_file_param{prj__chemical_system__kinetic_reactants}
+ config.getConfigSubtreeOptional("kinetic_reactants"), mesh,
+ chemical_system_map);
+
+ // equilibrium reactants
+ auto equilibrium_reactants = createEquilibriumReactants(
+ //! \ogs_file_param{prj__chemical_system__equilibrium_reactants}
+ config.getConfigSubtreeOptional("equilibrium_reactants"), mesh,
+ chemical_system_map);
+
+ return std::make_unique<ChemicalSystem>(std::move(aqueous_solution),
+ std::move(kinetic_reactants),
+ std::move(equilibrium_reactants));
+}
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.h b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..59962d127534d1407aaa6c0b6c5eca093220276b
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/CreateChemicalSystem.h
@@ -0,0 +1,38 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+} // namespace MeshLib
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+struct ChemicalSystem;
+
+std::unique_ptr<ChemicalSystem> createChemicalSystem(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh,
+ MeshLib::PropertyVector<std::size_t> const& chemical_system_map);
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateOutput.cpp b/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
index b34b7441de62e752121ca7a35ce885069992fe2a..70561477c1aa51b16fc56697e3e376876b9f6cba 100644
--- a/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
@@ -12,8 +12,8 @@
#include <numeric>
#include "AqueousSolution.h"
+#include "ChemicalSystem.h"
#include "CreateOutput.h"
-#include "EquilibriumReactant.h"
#include "KineticReactant.h"
#include "UserPunch.h"
@@ -22,13 +22,14 @@ namespace ChemistryLib
namespace PhreeqcIOData
{
std::unique_ptr<Output> createOutput(
- std::vector<Component> const& components,
- std::vector<EquilibriumReactant> const& equilibrium_reactants,
- std::vector<KineticReactant> const& kinetic_reactants,
+ ChemicalSystem const& chemical_system,
std::unique_ptr<UserPunch> const& user_punch,
bool const use_high_precision,
std::string const& project_file_name)
{
+ auto const& components = chemical_system.aqueous_solution->components;
+ auto const& equilibrium_reactants = chemical_system.equilibrium_reactants;
+ auto const& kinetic_reactants = chemical_system.kinetic_reactants;
// Mark which phreeqc output items will be held.
std::vector<OutputItem> accepted_items{{"pH", ItemType::pH},
{"pe", ItemType::pe}};
diff --git a/ChemistryLib/PhreeqcIOData/CreateOutput.h b/ChemistryLib/PhreeqcIOData/CreateOutput.h
index 1f542fc764dff989f1c2b555e6d0d4d3c0564c60..25663db4be0cbd814a41bb156d481453347f5c65 100644
--- a/ChemistryLib/PhreeqcIOData/CreateOutput.h
+++ b/ChemistryLib/PhreeqcIOData/CreateOutput.h
@@ -18,15 +18,11 @@ namespace ChemistryLib
namespace PhreeqcIOData
{
struct Output;
-struct Component;
-struct EquilibriumReactant;
-struct KineticReactant;
+struct ChemicalSystem;
struct UserPunch;
std::unique_ptr<Output> createOutput(
- std::vector<Component> const& components,
- std::vector<EquilibriumReactant> const& equilibrium_reactants,
- std::vector<KineticReactant> const& kinetic_reactants,
+ ChemicalSystem const& chemical_system,
std::unique_ptr<UserPunch> const& user_punch,
bool const use_high_precision,
std::string const& project_file_name);
diff --git a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp
index 44c0fce8422b462bcd12a88f9489c4d39f877ce5..3d9d25577b47325b589bfc6f812ffa7d22a70b1f 100644
--- a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.cpp
@@ -17,7 +17,7 @@ namespace ChemistryLib
namespace PhreeqcIOData
{
std::vector<Component> createSolutionComponents(
- BaseLib::ConfigTree const& config)
+ BaseLib::ConfigTree const& config, std::size_t const num_chemical_systems)
{
std::vector<Component> components;
//! \ogs_file_param{prj__chemical_system__solution__components}
@@ -27,7 +27,7 @@ std::vector<Component> createSolutionComponents(
//! \ogs_file_param{prj__chemical_system__solution__components__component}
comp_config.getConfigParameterList<std::string>("component"))
{
- components.emplace_back(component_name);
+ components.emplace_back(component_name, num_chemical_systems);
}
return components;
diff --git a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h
index 417450e9dac2319acd7b65572483082a49ecefb1..d745e90979c74df7e02274a90684b61ba6d01ddc 100644
--- a/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h
+++ b/ChemistryLib/PhreeqcIOData/CreateSolutionComponent.h
@@ -24,6 +24,6 @@ namespace PhreeqcIOData
struct Component;
std::vector<Component> createSolutionComponents(
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config, std::size_t const num_chemical_systems);
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/MathLib/LinAlg/MatrixVectorTraits.cpp b/MathLib/LinAlg/MatrixVectorTraits.cpp
index cac6b4f51f18ee50b03b0831c7a9b5f91eff64e7..ffc7174c86795620aba0c56e4dcba4538b3c6f84 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.cpp
+++ b/MathLib/LinAlg/MatrixVectorTraits.cpp
@@ -82,6 +82,13 @@ newInstance(MatrixSpecifications const& spec)
}
}
+std::unique_ptr<PETScVector> MatrixVectorTraits<PETScVector>::newInstance(
+ PETScVector::IndexType const length)
+{
+ auto const is_global_size = true;
+
+ return std::make_unique<PETScVector>(length, is_global_size);
+}
} // namespace MathLib
@@ -139,6 +146,11 @@ newInstance(MatrixSpecifications const& spec)
return std::make_unique<EigenVector>(spec.nrows);
}
+std::unique_ptr<EigenVector> MatrixVectorTraits<EigenVector>::newInstance(
+ Eigen::SparseMatrix<double>::Index const length)
+{
+ return std::make_unique<EigenVector>(length);
+}
} // namespace MathLib
#endif // defined(OGS_USE_EIGEN)
diff --git a/MathLib/LinAlg/MatrixVectorTraits.h b/MathLib/LinAlg/MatrixVectorTraits.h
index 608dce16719c087e939468866faafeb978e1c224..4e6f2737b5a01f8e68ad2371f8d16c6b8e6433ec 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.h
+++ b/MathLib/LinAlg/MatrixVectorTraits.h
@@ -19,15 +19,18 @@ template<typename Matrix>
struct MatrixVectorTraits;
}
-#define SPECIALIZE_MATRIX_VECTOR_TRAITS(MATVEC, IDX) \
- template<> struct MatrixVectorTraits<MATVEC> { \
- using Index = IDX; \
- static std::unique_ptr<MATVEC> newInstance(); \
- static std::unique_ptr<MATVEC> newInstance(MATVEC const& A); \
- static std::unique_ptr<MATVEC> newInstance(MatrixSpecifications const& spec); \
+#define SPECIALIZE_MATRIX_VECTOR_TRAITS(MATVEC, IDX) \
+ template <> \
+ struct MatrixVectorTraits<MATVEC> \
+ { \
+ using Index = IDX; \
+ static std::unique_ptr<MATVEC> newInstance(); \
+ static std::unique_ptr<MATVEC> newInstance(MATVEC const& A); \
+ static std::unique_ptr<MATVEC> newInstance( \
+ MatrixSpecifications const& spec); \
+ static std::unique_ptr<MATVEC> newInstance(Index const length); \
};
-
#ifdef USE_PETSC
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"