15 #include <boost/algorithm/string.hpp>
41 namespace PhreeqcIOData
45 template <
typename DataBlock>
47 std::vector<DataBlock>
const& data_blocks)
49 std::copy(data_blocks.begin(), data_blocks.end(),
50 std::ostream_iterator<DataBlock>(os));
59 auto& components = aqueous_solution.
components;
60 for (
unsigned component_id = 0; component_id < components.size();
63 components[component_id].amount->set(chemical_system_id,
64 concentrations[component_id]);
68 aqueous_solution.
pH->set(chemical_system_id, concentrations.back());
71 template <
typename Reactant>
80 auto const& solid_constituent = solid_phase.
component(reactant.name);
82 if (solid_constituent.hasProperty(
87 .template initialValue<double>(pos, t);
89 (*reactant.molality)[chemical_system_id] =
molality;
90 (*reactant.molality_prev)[chemical_system_id] =
molality;
97 .template initialValue<double>(pos, t);
101 (*reactant.volume_fraction_prev)[chemical_system_id] =
volume_fraction;
105 .template initialValue<double>(pos, t);
107 auto const porosity =
109 .template initialValue<double>(pos, t);
113 .template initialValue<double>(pos, t);
115 (*reactant.molality)[chemical_system_id] =
118 (*reactant.molality_prev)[chemical_system_id] =
119 (*reactant.molality)[chemical_system_id];
123 template <
typename Reactant>
130 double const t,
double const dt)
132 auto const& solid_constituent = solid_phase.
component(reactant.name);
134 if (solid_constituent.hasProperty(
137 (*reactant.molality_prev)[chemical_system_id] =
138 (*reactant.molality)[chemical_system_id];
144 (*reactant.volume_fraction)[chemical_system_id];
146 (*reactant.volume_fraction_prev)[chemical_system_id] =
147 (*reactant.volume_fraction)[chemical_system_id];
151 .template value<double>(vars, pos, t, dt);
153 auto const porosity = std::get<double>(
158 .template value<double>(vars, pos, t, dt);
160 (*reactant.molality)[chemical_system_id] =
163 (*reactant.molality_prev)[chemical_system_id] =
164 (*reactant.molality)[chemical_system_id];
167 template <
typename Exchanger>
174 auto const& solid_constituent = solid_phase.
component(exchanger.name);
178 .template initialValue<double>(pos, t);
180 (*exchanger.molality)[chemical_system_id] =
molality;
183 template <
typename Reactant>
188 double const porosity,
double const t,
191 auto const& solid_phase = medium.
phase(
"Solid");
192 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
196 auto const liquid_density =
198 .template value<double>(vars, pos, t, dt);
200 auto const& solid_constituent = solid_phase.component(reactant.name);
202 if (solid_constituent.hasProperty(
210 .template value<double>(vars, pos, t, dt);
212 (*reactant.volume_fraction)[chemical_system_id] +=
213 ((*reactant.molality)[chemical_system_id] -
214 (*reactant.molality_prev)[chemical_system_id]) *
218 template <
typename Reactant>
224 auto const& solid_phase = medium.
phase(
"Solid");
226 auto const& solid_constituent = solid_phase.
component(reactant.name);
228 if (solid_constituent.hasProperty(
234 porosity -= ((*reactant.volume_fraction)[chemical_system_id] -
235 (*reactant.volume_fraction_prev)[chemical_system_id]);
238 template <
typename Reactant>
240 Reactant
const& reactant,
241 std::vector<GlobalIndexType>
const& chemical_system_indices)
243 double const sum = std::accumulate(
244 chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
246 { return s + (*reactant.molality)[id]; });
247 return sum / chemical_system_indices.size();
252 std::string
const& project_file_name,
253 std::string&& database,
254 std::unique_ptr<ChemicalSystem>&& chemical_system,
255 std::vector<ReactionRate>&& reaction_rates,
256 std::vector<SurfaceSite>&& surface,
257 std::unique_ptr<UserPunch>&& user_punch,
258 std::unique_ptr<Output>&& output,
259 std::unique_ptr<Dump>&& dump,
262 _phreeqc_input_file(project_file_name +
"_phreeqc.inp"),
263 _database(std::move(database)),
264 _chemical_system(std::move(chemical_system)),
265 _reaction_rates(std::move(reaction_rates)),
266 _surface(std::move(surface)),
267 _user_punch(std::move(user_punch)),
268 _output(std::move(output)),
269 _dump(std::move(dump)),
270 _knobs(std::move(knobs))
276 "Failed to initialize phreeqc instance, due to lack of memory.");
283 "Failed in loading the specified thermodynamic database file: "
291 "Failed to fly the flag for the specified file {:s} where phreeqc "
292 "will write output.",
293 _output->basic_output_setups.output_file);
318 std::vector<double>
const& concentrations,
327 auto const& solid_phase = medium.
phase(
"Solid");
328 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
333 solid_phase, liquid_phase, medium, pos, t);
339 solid_phase, liquid_phase, medium, pos, t);
350 std::vector<double>
const& concentrations,
359 auto const& solid_phase = medium->
phase(
"Solid");
360 auto const& liquid_phase = medium->
phase(
"AqueousLiquid");
365 liquid_phase, vars, pos, t, dt);
371 solid_phase, liquid_phase, vars, pos, t, dt);
387 std::vector<GlobalVector*> int_pt_process_solutions;
388 int_pt_process_solutions.reserve(aqueous_solution->components.size() + 1);
390 std::transform(aqueous_solution->components.begin(),
391 aqueous_solution->components.end(),
392 std::back_inserter(int_pt_process_solutions),
393 [](
auto const&
c) { return c.amount.get(); });
395 int_pt_process_solutions.push_back(aqueous_solution->pH.get());
397 return int_pt_process_solutions;
401 int const component_id,
GlobalIndexType const chemical_system_id)
const
404 auto& components = aqueous_solution.components;
405 auto const&
pH = *aqueous_solution.pH;
407 return component_id < static_cast<int>(components.size())
408 ? components[component_id].amount->get(chemical_system_id)
409 :
pH.get(chemical_system_id);
419 auto const& dump_file =
_dump->dump_file;
420 std::ifstream in(dump_file);
423 OGS_FATAL(
"Could not open phreeqc dump file '{:s}'.", dump_file);
430 OGS_FATAL(
"Error when reading phreeqc dump file '{:s}'", dump_file);
443 OGS_FATAL(
"Could not open file '{:s}' for writing phreeqc inputs.",
447 out << std::scientific
448 << std::setprecision(std::numeric_limits<double>::digits10);
449 out << (*
this << dt);
453 OGS_FATAL(
"Failed in generating phreeqc input file '{:s}'.",
462 os << phreeqc_io.
_knobs <<
"\n";
464 os << *phreeqc_io.
_output <<
"\n";
469 os << *user_punch <<
"\n";
473 if (!reaction_rates.empty())
477 os << reaction_rates <<
"\n";
480 for (std::size_t chemical_system_id = 0;
482 ++chemical_system_id)
484 os <<
"SOLUTION " << chemical_system_id + 1 <<
"\n";
486 os, chemical_system_id);
488 auto const& dump = phreeqc_io.
_dump;
491 auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
492 if (!aqueous_solutions_prev.empty())
494 os << aqueous_solutions_prev[chemical_system_id] <<
"\n\n";
498 os <<
"USE solution none"
503 os <<
"USE solution " << chemical_system_id + 1 <<
"\n\n";
505 auto const& equilibrium_reactants =
507 if (!equilibrium_reactants.empty())
509 os <<
"EQUILIBRIUM_PHASES " << chemical_system_id + 1 <<
"\n";
510 for (
auto const& equilibrium_reactant : equilibrium_reactants)
512 equilibrium_reactant.print(os, chemical_system_id);
517 auto const& kinetic_reactants =
519 if (!kinetic_reactants.empty())
521 os <<
"KINETICS " << chemical_system_id + 1 <<
"\n";
522 for (
auto const& kinetic_reactant : kinetic_reactants)
524 kinetic_reactant.print(os, chemical_system_id);
526 os <<
"-steps " << phreeqc_io.
_dt <<
"\n"
530 auto const& surface = phreeqc_io.
_surface;
531 if (!surface.empty())
533 os <<
"SURFACE " << chemical_system_id + 1 <<
"\n";
534 std::size_t aqueous_solution_id =
535 dump->aqueous_solutions_prev.empty()
536 ? chemical_system_id + 1
538 os <<
"-equilibrate with solution " << aqueous_solution_id <<
"\n";
539 os <<
"-sites_units DENSITY"
541 os << surface <<
"\n";
542 os <<
"SAVE solution " << chemical_system_id + 1 <<
"\n";
546 if (!exchangers.empty())
548 os <<
"EXCHANGE " << chemical_system_id + 1 <<
"\n";
549 std::size_t
const aqueous_solution_id =
550 dump->aqueous_solutions_prev.empty()
551 ? chemical_system_id + 1
553 os <<
"-equilibrate with solution " << aqueous_solution_id <<
"\n";
554 for (
auto const& exchanger : exchangers)
556 exchanger.print(os, chemical_system_id);
558 os <<
"SAVE solution " << chemical_system_id + 1 <<
"\n";
565 auto const& dump = phreeqc_io.
_dump;
576 INFO(
"Phreeqc: Executing chemical calculation.");
581 "Failed in performing speciation calculation with the generated "
582 "phreeqc input file '{:s}'.",
589 auto const& basic_output_setups =
_output->basic_output_setups;
590 auto const& phreeqc_result_file = basic_output_setups.output_file;
591 DBUG(
"Reading phreeqc results from file '{:s}'.", phreeqc_result_file);
592 std::ifstream in(phreeqc_result_file);
596 OGS_FATAL(
"Could not open phreeqc result file '{:s}'.",
597 phreeqc_result_file);
604 OGS_FATAL(
"Error when reading phreeqc result file '{:s}'",
605 phreeqc_result_file);
614 in.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
617 auto const& output = *phreeqc_io.
_output;
618 auto const& dropped_item_ids = output.dropped_item_ids;
620 auto const& surface = phreeqc_io.
_surface;
623 if (!surface.empty() && !exchangers.empty())
626 "Using surface and exchange reactions simultaneously is not "
627 "supported at the moment");
630 int const num_skipped_lines = surface.empty() && exchangers.empty() ? 1 : 2;
632 auto& equilibrium_reactants =
636 for (std::size_t chemical_system_id = 0;
638 ++chemical_system_id)
641 for (
int i = 0; i < num_skipped_lines; ++i)
643 in.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
647 if (!std::getline(in, line))
650 "Error when reading calculation result of Solution {:d} "
651 "after the reaction.",
655 std::vector<double> accepted_items;
656 std::vector<std::string> items;
657 boost::trim_if(line, boost::is_any_of(
"\t "));
658 boost::algorithm::split(items, line, boost::is_any_of(
"\t "),
659 boost::token_compress_on);
660 for (
int item_id = 0; item_id < static_cast<int>(items.size());
663 if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
664 item_id) == dropped_item_ids.end())
669 value = std::stod(items[item_id]);
671 catch (
const std::invalid_argument& e)
674 "Invalid argument. Could not convert string '{:s}' to "
675 "double for chemical system {:d}, column {:d}. "
676 "Exception '{:s}' was thrown.",
677 items[item_id], chemical_system_id + 1, item_id,
680 catch (
const std::out_of_range& e)
683 "Out of range error. Could not convert string "
684 "'{:s}' to double for chemical system {:d}, column "
685 "{:d}. Exception '{:s}' was thrown.",
686 items[item_id], chemical_system_id + 1, item_id,
689 accepted_items.push_back(value);
692 assert(accepted_items.size() == output.accepted_items.size());
695 auto& components = aqueous_solution->components;
697 for (
int item_id = 0; item_id < static_cast<int>(accepted_items.size());
700 auto const& accepted_item = output.accepted_items[item_id];
701 auto const& item_name = accepted_item.name;
703 auto compare_by_name = [&item_name](
auto const& item)
704 {
return item.name == item_name; };
706 switch (accepted_item.item_type)
711 aqueous_solution->pH->set(
713 std::pow(10, -accepted_items[item_id]));
719 (*aqueous_solution->pe)[chemical_system_id] =
720 accepted_items[item_id];
727 components.begin(), components.end(), compare_by_name,
728 "Could not find component '" + item_name +
"'.");
729 component.amount->set(chemical_system_id,
730 accepted_items[item_id]);
737 equilibrium_reactants.begin(),
738 equilibrium_reactants.end(), compare_by_name,
739 "Could not find equilibrium reactant '" + item_name +
741 (*equilibrium_reactant.molality)[chemical_system_id] =
742 accepted_items[item_id];
749 kinetic_reactants.begin(), kinetic_reactants.end(),
751 "Could not find kinetic reactant '" + item_name +
"'.");
752 (*kinetic_reactant.molality)[chemical_system_id] =
753 accepted_items[item_id];
759 auto& secondary_variables = user_punch->secondary_variables;
762 secondary_variables.begin(), secondary_variables.end(),
764 "Could not find secondary variable '" + item_name +
766 (*secondary_variable.value)[chemical_system_id] =
767 accepted_items[item_id];
779 std::vector<std::string> component_names;
781 std::transform(components.begin(), components.end(),
782 std::back_inserter(component_names),
783 [](
auto const&
c) { return c.name; });
785 component_names.push_back(
"H");
787 return component_names;
794 double const t,
double const dt)
799 medium, pos, porosity, t, dt);
805 medium, pos, porosity, t, dt);
828 std::size_t
const ele_id,
829 std::vector<GlobalIndexType>
const& chemical_system_indices)
833 (*kinetic_reactant.mesh_prop_molality)[ele_id] =
839 (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
841 chemical_system_indices);
GlobalMatrix::IndexType GlobalIndexType
void INFO(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
Definition of the Mesh class.
std::vector< GlobalIndexType > chemical_system_index_map
std::vector< SurfaceSite > const _surface
std::string const _phreeqc_input_file
PhreeqcIO(GlobalLinearSolver &linear_solver, std::string const &project_file_name, std::string &&database, std::unique_ptr< ChemicalSystem > &&chemical_system, std::vector< ReactionRate > &&reaction_rates, std::vector< SurfaceSite > &&surface, std::unique_ptr< UserPunch > &&user_punch, std::unique_ptr< Output > &&output, std::unique_ptr< Dump > &&dump, Knobs &&knobs)
double getConcentration(int const component_id, GlobalIndexType const chemical_system_id) const override
void setAqueousSolutionsPrevFromDumpFile() override
std::unique_ptr< ChemicalSystem > _chemical_system
std::vector< ReactionRate > const _reaction_rates
std::unique_ptr< UserPunch > _user_punch
const int phreeqc_instance_id
void writeInputsToFile(double const dt)
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
std::unique_ptr< Dump > const _dump
void initialize() override
std::unique_ptr< Output > const _output
void readOutputsFromFile()
std::vector< GlobalVector * > getIntPtProcessSolutions() const override
std::size_t _num_chemical_systems
void computeSecondaryVariable(std::size_t const ele_id, std::vector< GlobalIndexType > const &chemical_system_indices) 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 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 executeSpeciationCalculation(double const dt) override
std::string const _database
std::vector< std::string > const getComponentList() const override
void updatePorosityPostReaction(GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, double &porosity) override
Phase const & phase(std::size_t index) const
Component const & component(std::size_t const &index) const
std::iterator_traits< InputIt >::reference findElementOrError(InputIt begin, InputIt end, Predicate predicate, std::string const &error="")
void initializeReactantMolality(Reactant &reactant, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Phase const &solid_phase, MaterialPropertyLib::Phase const &liquid_phase, MaterialPropertyLib::Medium const &medium, ParameterLib::SpatialPosition const &pos, double const t)
void initializeExchangerMolality(Exchanger &exchanger, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Phase const &solid_phase, ParameterLib::SpatialPosition const &pos, double const t)
void setPorosityPostReaction(Reactant &reactant, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, double &porosity)
void setReactantMolality(Reactant &reactant, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Phase const &solid_phase, MaterialPropertyLib::Phase const &liquid_phase, MaterialPropertyLib::VariableArray const &vars, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
void setAqueousSolution(std::vector< double > const &concentrations, GlobalIndexType const &chemical_system_id, AqueousSolution &aqueous_solution)
static double averageReactantMolality(Reactant const &reactant, std::vector< GlobalIndexType > const &chemical_system_indices)
void updateReactantVolumeFraction(Reactant &reactant, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, ParameterLib::SpatialPosition const &pos, double const porosity, double const t, double const dt)
std::istream & operator>>(std::istream &in, PhreeqcIO &phreeqc_io)
std::ostream & operator<<(std::ostream &os, PhreeqcIO const &phreeqc_io)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
void copy(PETScVector const &x, PETScVector &y)
double fluid_density(const double p, const double T, const double x)
std::vector< Component > components
std::unique_ptr< GlobalVector > pH