diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index d1d7683d4fa965789946c9add884499d7cd295e2..3f4a3d28253d3d816437f588f6d6f842b397435c 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -17,9 +17,11 @@
#include "PhreeqcIOData/CreateEquilibriumPhase.h"
#include "PhreeqcIOData/CreateKineticReactant.h"
#include "PhreeqcIOData/CreateOutput.h"
+#include "PhreeqcIOData/CreateUserPunch.h"
#include "PhreeqcIOData/EquilibriumPhase.h"
#include "PhreeqcIOData/KineticReactant.h"
#include "PhreeqcIOData/ReactionRate.h"
+#include "PhreeqcIOData/UserPunch.h"
#include "PhreeqcKernel.h"
#include "PhreeqcKernelData/AqueousSolution.h"
#include "PhreeqcKernelData/CreateAqueousSolution.h"
@@ -93,6 +95,11 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
auto const project_file_name = BaseLib::joinPaths(
output_directory,
BaseLib::extractBaseNameWithoutExtension(config.getProjectFileName()));
+
+ // user punch
+ auto user_punch = PhreeqcIOData::createUserPunch(
+ //! \ogs_file_param{prj__chemical_system__user_punch}
+ config.getConfigSubtreeOptional("user_punch"), mesh);
auto output = PhreeqcIOData::createOutput(
components, equilibrium_phases, kinetic_reactants, project_file_name);
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 5fd72adc44353fad7041c073254d03e4e6850e67..60c5f0c49931ed7bc7cd659cf6342d22c37d01f5 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -370,6 +370,8 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
auto& components = aqueous_solution.components;
auto& equilibrium_phases = phreeqc_io._equilibrium_phases;
auto& kinetic_reactants = phreeqc_io._kinetic_reactants;
+ auto& user_punch = phreeqc_io._user_punch;
+ auto& secondary_variables = user_punch->secondary_variables;
for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
++item_id)
{
@@ -426,6 +428,18 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
accepted_items[item_id];
break;
}
+ case ItemType::SecondaryVariable:
+ {
+ // Update values of secondary variables
+ auto& secondary_variable = BaseLib::findElementOrError(
+ secondary_variables.begin(), secondary_variables.end(),
+ compare_by_name,
+ "Could not find secondary variable'" + item_name +
+ "'.");
+ (*secondary_variable.value)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
}
}
}
diff --git a/ChemistryLib/PhreeqcIOData/CreateOutput.cpp b/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
index 7f2d5f9b7457e00673decdcff4d47819166779b5..669967e7b739f1bce793df45518c92baa736c299 100644
--- a/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateOutput.cpp
@@ -15,6 +15,7 @@
#include "CreateOutput.h"
#include "EquilibriumPhase.h"
#include "KineticReactant.h"
+#include "UserPunch.h"
namespace ChemistryLib
{
@@ -24,6 +25,7 @@ std::unique_ptr<Output> createOutput(
std::vector<Component> const& components,
std::vector<EquilibriumPhase> const& equilibrium_phases,
std::vector<KineticReactant> const& kinetic_reactants,
+ std::unique_ptr<UserPunch> const& user_punch,
std::string const& project_file_name)
{
// Mark which phreeqc output items will be held.
@@ -48,6 +50,13 @@ std::unique_ptr<Output> createOutput(
kinetic_reactant.item_type);
}
+ if (user_punch)
+ {
+ auto const& secondary_variables = user_punch->secondary_variables;
+ std::transform(secondary_variables.begin(), secondary_variables.end(),
+ std::back_inserter(accepted_items), accepted_item);
+ }
+
// Record ids of which phreeqc output items will be dropped.
BasicOutputSetups basic_output_setups(project_file_name);
auto const num_dropped_basic_items =
diff --git a/ChemistryLib/PhreeqcIOData/CreateOutput.h b/ChemistryLib/PhreeqcIOData/CreateOutput.h
index 3ef3993bb1a8384803567e0cfb708d02e1c56e0f..ce90837555837e1dbd346e8713e66a0ac876aaec 100644
--- a/ChemistryLib/PhreeqcIOData/CreateOutput.h
+++ b/ChemistryLib/PhreeqcIOData/CreateOutput.h
@@ -21,11 +21,13 @@ struct Output;
struct Component;
struct EquilibriumPhase;
struct KineticReactant;
+struct UserPunch;
std::unique_ptr<Output> createOutput(
std::vector<Component> const& components,
std::vector<EquilibriumPhase> const& equilibrium_phases,
std::vector<KineticReactant> const& kinetic_reactants,
+ std::unique_ptr<UserPunch> const& user_punch,
std::string const& project_file_name);
} // namespace PhreeqcIOData
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateUserPunch.cpp b/ChemistryLib/PhreeqcIOData/CreateUserPunch.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9fd59e131dacca96b98ca4bf5307e58aba2308d4
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/CreateUserPunch.cpp
@@ -0,0 +1,56 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <boost/optional/optional.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
+#include "UserPunch.h"
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+std::unique_ptr<UserPunch> createUserPunch(
+ boost::optional<BaseLib::ConfigTree> const& config,
+ MeshLib::Mesh const& mesh)
+{
+ if (!config)
+ return nullptr;
+
+ std::vector<SecondaryVariable> secondary_variables;
+ for (auto const& variable_name :
+ //! \ogs_file_param{prj__chemical_system__user_punch__headline}
+ config->getConfigParameter<std::vector<std::string>>("headline"))
+ {
+ auto value = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ variable_name,
+ MeshLib::MeshItemType::Node,
+ 1);
+ std::fill(std::begin(*value),
+ std::end(*value),
+ std::numeric_limits<double>::quiet_NaN());
+
+ secondary_variables.emplace_back(variable_name, value);
+ }
+
+ std::vector<std::string> statements;
+ for (auto const& statement :
+ //! \ogs_file_param{prj__chemical_system__user_punch__statement}
+ config->getConfigParameterList<std::string>("statement"))
+ {
+ statements.emplace_back(statement);
+ }
+
+ return std::make_unique<UserPunch>(std::move(secondary_variables),
+ std::move(statements));
+}
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateUserPunch.h b/ChemistryLib/PhreeqcIOData/CreateUserPunch.h
new file mode 100644
index 0000000000000000000000000000000000000000..48f3f7762b8a6438bbdfffb336ad88fc5f82220e
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/CreateUserPunch.h
@@ -0,0 +1,35 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <boost/optional/optional_fwd.hpp>
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+struct UserPunch;
+
+std::unique_ptr<UserPunch> createUserPunch(
+ boost::optional<BaseLib::ConfigTree> const& config,
+ MeshLib::Mesh const& mesh);
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/Output.h b/ChemistryLib/PhreeqcIOData/Output.h
index 69d1c808f0600d4a9112bd9ea5590ef8af6743f4..8b1322ecbde541803922ed55876dec8d32cf54d3 100644
--- a/ChemistryLib/PhreeqcIOData/Output.h
+++ b/ChemistryLib/PhreeqcIOData/Output.h
@@ -63,7 +63,8 @@ enum class ItemType
pe,
Component,
EquilibriumPhase,
- KineticReactant
+ KineticReactant,
+ SecondaryVariable
};
struct OutputItem
diff --git a/ChemistryLib/PhreeqcIOData/UserPunch.cpp b/ChemistryLib/PhreeqcIOData/UserPunch.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..07a4c736ba1b97b0748af9f789f735c3bff2f4d9
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/UserPunch.cpp
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "UserPunch.h"
+#include "BaseLib/ConfigTreeUtil.h"
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+std::ostream& operator<<(std::ostream& os, UserPunch const& user_punch)
+{
+ os << "USER_PUNCH" << "\n";
+ os << "-headings ";
+ auto const& secondary_variables = user_punch.secondary_variables;
+ for (auto& secondary_variable : secondary_variables)
+ os << secondary_variable.name.c_str() << " ";
+ os << "\n";
+
+ os << "-start" << "\n";
+ int line_number = 1;
+ for (auto const& statement : user_punch.statements)
+ {
+ os << line_number << " " << statement.c_str() << "\n";
+ ++line_number;
+ }
+ os << "-end" << "\n";
+
+ return os;
+}
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/UserPunch.h b/ChemistryLib/PhreeqcIOData/UserPunch.h
new file mode 100644
index 0000000000000000000000000000000000000000..5a1ceb425092a1e0cd4a2b0728c77fa4154dbf6a
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/UserPunch.h
@@ -0,0 +1,56 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <iosfwd>
+#include <string>
+#include <vector>
+
+#include "MeshLib/PropertyVector.h"
+#include "Output.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+struct SecondaryVariable
+{
+ SecondaryVariable(std::string const name_,
+ MeshLib::PropertyVector<double>* value_)
+ : name(name_), value(value_)
+ {
+ }
+
+ std::string const name;
+ MeshLib::PropertyVector<double>* value;
+ static const ItemType item_type = ItemType::SecondaryVariable;
+};
+
+struct UserPunch
+{
+ UserPunch(std::vector<SecondaryVariable>&& secondary_variables_,
+ std::vector<std::string>&& statements_)
+ : secondary_variables(secondary_variables_), statements(statements_)
+ {
+ }
+
+ friend std::ostream& operator<<(std::ostream& os,
+ UserPunch const& user_punch);
+
+ std::vector<SecondaryVariable> secondary_variables;
+ std::vector<std::string> const statements;
+};
+} // namespace PhreeqcIOData
+} // namespace ChemistryLib