OGS
PhreeqcIO.cpp
Go to the documentation of this file.
1
11#include "PhreeqcIO.h"
12
13#include <IPhreeqc.h>
14
15#include <boost/algorithm/string.hpp>
16#include <cmath>
17#include <fstream>
18#include <iomanip>
19#include <numeric>
20
21#include "BaseLib/Algorithm.h"
26#include "MeshLib/Mesh.h"
29#include "PhreeqcIOData/Dump.h"
33#include "PhreeqcIOData/Knobs.h"
38
39namespace ChemistryLib
40{
41namespace PhreeqcIOData
42{
43namespace
44{
45template <class... Ts>
46struct overloaded : Ts...
47{
48 using Ts::operator()...;
49};
50template <class... Ts>
51overloaded(Ts...) -> overloaded<Ts...>;
52
53template <typename DataBlock>
54std::ostream& operator<<(std::ostream& os,
55 std::vector<DataBlock> const& data_blocks)
56{
57 std::copy(data_blocks.begin(), data_blocks.end(),
58 std::ostream_iterator<DataBlock>(os));
59 return os;
60}
61
62void setAqueousSolution(std::vector<double> const& concentrations,
63 GlobalIndexType const& chemical_system_id,
64 AqueousSolution& aqueous_solution)
65{
66 GlobalIndexType const offset = aqueous_solution.pH->getRangeBegin();
67 GlobalIndexType const global_index = offset + chemical_system_id;
68
69 // components
70 auto& components = aqueous_solution.components;
71 for (unsigned component_id = 0; component_id < components.size();
72 ++component_id)
73 {
75 *components[component_id].amount);
76 components[component_id].amount->set(global_index,
77 concentrations[component_id]);
78 }
79
80 // pH
82 aqueous_solution.pH->set(global_index, concentrations.back());
83}
84
85template <typename Reactant>
86void initializeReactantMolality(Reactant& reactant,
87 GlobalIndexType const& chemical_system_id,
88 MaterialPropertyLib::Phase const& solid_phase,
89 MaterialPropertyLib::Phase const& liquid_phase,
90 MaterialPropertyLib::Medium const& medium,
92 double const t)
93{
94 auto const& solid_constituent = solid_phase.component(reactant.name);
95
96 if (solid_constituent.hasProperty(
98 {
99 auto const molality =
101 .template initialValue<double>(pos, t);
102
103 (*reactant.molality)[chemical_system_id] = molality;
104 (*reactant.molality_prev)[chemical_system_id] = molality;
105 }
106 else
107 {
108 auto const volume_fraction =
109 solid_constituent
111 .template initialValue<double>(pos, t);
112
113 (*reactant.volume_fraction)[chemical_system_id] = volume_fraction;
114
115 (*reactant.volume_fraction_prev)[chemical_system_id] = volume_fraction;
116
117 auto const fluid_density =
119 .template initialValue<double>(pos, t);
120
121 auto const porosity =
123 .template initialValue<double>(pos, t);
124
125 auto const molar_volume =
127 .template initialValue<double>(pos, t);
128
129 (*reactant.molality)[chemical_system_id] =
130 volume_fraction / fluid_density / porosity / molar_volume;
131
132 (*reactant.molality_prev)[chemical_system_id] =
133 (*reactant.molality)[chemical_system_id];
134 }
135}
136
137template <typename Reactant>
138void setReactantMolality(Reactant& reactant,
139 GlobalIndexType const& chemical_system_id,
140 MaterialPropertyLib::Phase const& solid_phase,
141 MaterialPropertyLib::Phase const& liquid_phase,
144 double const t, double const dt)
145{
146 auto const& solid_constituent = solid_phase.component(reactant.name);
147
148 if (solid_constituent.hasProperty(
150 {
151 (*reactant.molality_prev)[chemical_system_id] =
152 (*reactant.molality)[chemical_system_id];
153
154 return;
155 }
156
157 auto const volume_fraction =
158 (*reactant.volume_fraction)[chemical_system_id];
159
160 (*reactant.volume_fraction_prev)[chemical_system_id] =
161 (*reactant.volume_fraction)[chemical_system_id];
162
163 auto const fluid_density =
165 .template value<double>(vars, pos, t, dt);
166
167 auto const molar_volume =
169 .template value<double>(vars, pos, t, dt);
170
171 (*reactant.molality)[chemical_system_id] =
172 volume_fraction / fluid_density / vars.porosity / molar_volume;
173
174 (*reactant.molality_prev)[chemical_system_id] =
175 (*reactant.molality)[chemical_system_id];
176}
177
178template <typename Site>
180 GlobalIndexType const& chemical_system_id,
181 MaterialPropertyLib::Phase const& solid_phase,
183 double const t)
184{
185 auto const& solid_constituent = solid_phase.component(site.name);
186
187 auto const molality =
189 .template initialValue<double>(pos, t);
190
191 (*site.molality)[chemical_system_id] = molality;
192}
193
194template <typename Reactant>
195void updateReactantVolumeFraction(Reactant& reactant,
196 GlobalIndexType const& chemical_system_id,
197 MaterialPropertyLib::Medium const& medium,
199 double const porosity, double const t,
200 double const dt)
201{
202 auto const& solid_phase = medium.phase("Solid");
203 auto const& liquid_phase = medium.phase("AqueousLiquid");
204
206
207 auto const liquid_density =
209 .template value<double>(vars, pos, t, dt);
210
211 auto const& solid_constituent = solid_phase.component(reactant.name);
212
213 if (solid_constituent.hasProperty(
215 {
216 return;
217 }
218
219 auto const molar_volume =
221 .template value<double>(vars, pos, t, dt);
222
223 (*reactant.volume_fraction)[chemical_system_id] +=
224 ((*reactant.molality)[chemical_system_id] -
225 (*reactant.molality_prev)[chemical_system_id]) *
226 liquid_density * porosity * molar_volume;
227}
228
229template <typename Reactant>
230void setPorosityPostReaction(Reactant& reactant,
231 GlobalIndexType const& chemical_system_id,
232 MaterialPropertyLib::Medium const& medium,
233 double& porosity)
234{
235 auto const& solid_phase = medium.phase("Solid");
236
237 auto const& solid_constituent = solid_phase.component(reactant.name);
238
239 if (solid_constituent.hasProperty(
241 {
242 return;
243 }
244
245 porosity -= ((*reactant.volume_fraction)[chemical_system_id] -
246 (*reactant.volume_fraction_prev)[chemical_system_id]);
247}
248
249template <typename Reactant>
251 Reactant const& reactant,
252 std::vector<GlobalIndexType> const& chemical_system_indices)
253{
254 double const sum = std::accumulate(
255 chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
256 [&](double const s, GlobalIndexType const id)
257 { return s + (*reactant.molality)[id]; });
258 return sum / chemical_system_indices.size();
259}
260} // namespace
261
262extern std::string specifyFileName(std::string const& project_file_name,
263 std::string const& file_extension);
264
266 GlobalLinearSolver& linear_solver,
267 std::string const& project_file_name,
268 std::string&& database,
269 std::unique_ptr<ChemicalSystem>&& chemical_system,
270 std::vector<ReactionRate>&& reaction_rates,
271 std::unique_ptr<UserPunch>&& user_punch,
272 std::unique_ptr<Output>&& output,
273 std::unique_ptr<Dump>&& dump,
274 Knobs&& knobs)
275 : ChemicalSolverInterface(mesh, linear_solver),
276 _phreeqc_input_file(specifyFileName(project_file_name, ".inp")),
277 _database(std::move(database)),
278 _chemical_system(std::move(chemical_system)),
279 _reaction_rates(std::move(reaction_rates)),
280 _user_punch(std::move(user_punch)),
281 _output(std::move(output)),
282 _dump(std::move(dump)),
283 _knobs(std::move(knobs))
284{
285 // initialize phreeqc instance
286 if (CreateIPhreeqc() != phreeqc_instance_id)
287 {
288 OGS_FATAL(
289 "Failed to initialize phreeqc instance, due to lack of memory.");
290 }
291
292 // load specified thermodynamic database
293 if (LoadDatabase(phreeqc_instance_id, _database.c_str()) != IPQ_OK)
294 {
295 OGS_FATAL(
296 "Failed in loading the specified thermodynamic database file: "
297 "{:s}.",
298 _database);
299 }
300
301 if (SetSelectedOutputFileOn(phreeqc_instance_id, 1) != IPQ_OK)
302 {
303 OGS_FATAL(
304 "Failed to fly the flag for the specified file {:s} where phreeqc "
305 "will write output.",
306 _output->basic_output_setups.output_file);
307 }
308
309 if (_dump)
310 {
311 // Chemical composition of the aqueous solution of last time step will
312 // be written into .dmp file once the second function argument is set to
313 // one.
314 SetDumpFileOn(phreeqc_instance_id, 1);
315 }
316}
317
319{
320 DestroyIPhreeqc(phreeqc_instance_id);
321}
322
324{
326
328
329 if (_user_punch)
330 {
332 }
333}
334
336 std::vector<double> const& concentrations,
337 GlobalIndexType const& chemical_system_id,
338 MaterialPropertyLib::Medium const& medium,
340 double const t)
341{
342 setAqueousSolution(concentrations, chemical_system_id,
343 *_chemical_system->aqueous_solution);
344
345 auto const& solid_phase = medium.phase("Solid");
346 auto const& liquid_phase = medium.phase("AqueousLiquid");
347
348 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
349 {
350 initializeReactantMolality(kinetic_reactant, chemical_system_id,
351 solid_phase, liquid_phase, medium, pos, t);
352 }
353
354 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
355 {
356 initializeReactantMolality(equilibrium_reactant, chemical_system_id,
357 solid_phase, liquid_phase, medium, pos, t);
358 }
359
360 for (auto& exchanger : _chemical_system->exchangers)
361 {
362 initializeSiteMolality(exchanger, chemical_system_id, solid_phase, pos,
363 t);
364 }
365
366 for (auto& surface_site : _chemical_system->surface)
367 {
368 if (auto const surface_site_ptr =
369 std::get_if<MoleBasedSurfaceSite>(&surface_site))
370 {
371 initializeSiteMolality(*surface_site_ptr, chemical_system_id,
372 solid_phase, pos, t);
373 }
374 }
375}
376
378 std::vector<double> const& concentrations,
379 GlobalIndexType const& chemical_system_id,
380 MaterialPropertyLib::Medium const* medium,
382 ParameterLib::SpatialPosition const& pos, double const t, double const dt)
383{
384 setAqueousSolution(concentrations, chemical_system_id,
385 *_chemical_system->aqueous_solution);
386
387 auto const& solid_phase = medium->phase("Solid");
388 auto const& liquid_phase = medium->phase("AqueousLiquid");
389
390 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
391 {
392 setReactantMolality(kinetic_reactant, chemical_system_id, solid_phase,
393 liquid_phase, vars, pos, t, dt);
394 }
395
396 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
397 {
398 setReactantMolality(equilibrium_reactant, chemical_system_id,
399 solid_phase, liquid_phase, vars, pos, t, dt);
400 }
401}
402
404{
406
407 callPhreeqc();
408
410}
411
413 int const component_id, GlobalIndexType const chemical_system_id) const
414{
415 auto const& aqueous_solution = *_chemical_system->aqueous_solution;
416 auto& components = aqueous_solution.components;
417 auto const& pH = *aqueous_solution.pH;
418
419 if (component_id < static_cast<int>(components.size()))
420 {
422 *components[component_id].amount);
423
424 return components[component_id].amount->get(chemical_system_id);
425 }
426
427 // pH
428 MathLib::LinAlg::setLocalAccessibleVector(*aqueous_solution.pH);
429
430 return pH.get(chemical_system_id);
431}
432
434{
435 if (!_dump)
436 {
437 return;
438 }
439
440 auto const& dump_file = _dump->dump_file;
441 std::ifstream in(dump_file);
442 if (!in)
443 {
444 // return if phreeqc dump file doesn't exist. This happens in
445 // the first time step when no dump file is provided by the user.
446 return;
447 }
448
449 _dump->readDumpFile(in, _num_chemical_systems);
450
451 if (!in)
452 {
453 OGS_FATAL("Error when reading phreeqc dump file '{:s}'", dump_file);
454 }
455
456 in.close();
457}
458
459void PhreeqcIO::writeInputsToFile(double const dt)
460{
461 DBUG("Writing phreeqc inputs into file '{:s}'.", _phreeqc_input_file);
462 std::ofstream out(_phreeqc_input_file, std::ofstream::out);
463
464 if (!out)
465 {
466 OGS_FATAL("Could not open file '{:s}' for writing phreeqc inputs.",
468 }
469
470 out << std::scientific
471 << std::setprecision(std::numeric_limits<double>::max_digits10);
472 out << (*this << dt);
473
474 if (!out)
475 {
476 OGS_FATAL("Failed in generating phreeqc input file '{:s}'.",
478 }
479
480 out.close();
481}
482
483std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
484{
485 bool const fixing_pe =
486 phreeqc_io._chemical_system->aqueous_solution->fixing_pe;
487 if (fixing_pe)
488 {
489 os << "PHASES\n"
490 << "Fix_pe\n"
491 << "e- = e-\n"
492 << "log_k 0.0\n\n";
493 }
494
495 os << phreeqc_io._knobs << "\n";
496
497 os << *phreeqc_io._output << "\n";
498
499 auto const& user_punch = phreeqc_io._user_punch;
500 if (user_punch)
501 {
502 os << *user_punch << "\n";
503 }
504
505 auto const& reaction_rates = phreeqc_io._reaction_rates;
506 if (!reaction_rates.empty())
507 {
508 os << "RATES\n";
509 os << reaction_rates << "\n";
510 }
511
512 for (std::size_t chemical_system_id = 0;
513 chemical_system_id < phreeqc_io._num_chemical_systems;
514 ++chemical_system_id)
515 {
516 os << "SOLUTION " << chemical_system_id + 1 << "\n";
517 phreeqc_io._chemical_system->aqueous_solution->print(
518 os, chemical_system_id);
519
520 auto const& dump = phreeqc_io._dump;
521 if (dump)
522 {
523 auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
524 // Print previous aqueous solution if it exists.
525 if (!aqueous_solutions_prev.empty())
526 {
527 os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
528 }
529 }
530
531 os << "USE solution none\n";
532 os << "END\n\n";
533
534 os << "USE solution " << chemical_system_id + 1 << "\n\n";
535
536 auto const& equilibrium_reactants =
537 phreeqc_io._chemical_system->equilibrium_reactants;
538 if (!equilibrium_reactants.empty() || fixing_pe)
539 {
540 os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
541 for (auto const& equilibrium_reactant : equilibrium_reactants)
542 {
543 equilibrium_reactant.print(os, chemical_system_id);
544 }
545 fixing_pe
546 ? os << "Fix_pe "
547 << -phreeqc_io._chemical_system->aqueous_solution->pe0
548 << " O2(g)\n\n"
549 : os << "\n";
550 }
551
552 auto const& kinetic_reactants =
553 phreeqc_io._chemical_system->kinetic_reactants;
554 if (!kinetic_reactants.empty())
555 {
556 os << "KINETICS " << chemical_system_id + 1 << "\n";
557 for (auto const& kinetic_reactant : kinetic_reactants)
558 {
559 kinetic_reactant.print(os, chemical_system_id);
560 }
561 os << "-steps " << phreeqc_io._dt << "\n\n";
562 }
563
564 auto const& surface = phreeqc_io._chemical_system->surface;
565 if (!surface.empty())
566 {
567 // To get the amount of surface species from the previous time step,
568 // an equilibration calculation with the previous aqueous solution
569 // is needed. The previous aqueous solution is saved using the
570 // PHREEQC keyword "DUMP" in the dump file and then stored as
571 // SOLUTION_RAW within aqueous_solutions_prev. Along with the
572 // PHREEQC keyword 'SURFACE', distinguish between the current and
573 // previous solutions, the previous solution number is added to the
574 // total number of chemical systems (_num_chemical_systems).
575 os << "SURFACE " << chemical_system_id + 1 << "\n";
576 std::size_t aqueous_solution_id =
577 dump->aqueous_solutions_prev.empty()
578 ? chemical_system_id + 1
579 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
580 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
581
582 // print unit
583 if (std::holds_alternative<DensityBasedSurfaceSite>(
584 surface.front()))
585 {
586 os << "-sites_units density\n";
587 }
588 else
589 {
590 os << "-sites_units absolute\n";
591 }
592
593 for (auto const& surface_site : surface)
594 {
595 std::visit(
596 overloaded{[&os](DensityBasedSurfaceSite const& s)
597 {
598 os << s.name << " " << s.site_density << " "
599 << s.specific_surface_area << " "
600 << s.mass << "\n";
601 },
602 [&os, chemical_system_id](
603 MoleBasedSurfaceSite const& s) {
604 os << s.name << " "
605 << (*s.molality)[chemical_system_id]
606 << "\n";
607 }},
608 surface_site);
609 }
610
611 // overlook the effect of the buildup of charges onto the surface
612 if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
613 {
614 os << "-no_edl\n";
615 }
616 // Here the PHREEQC "save" command is printed to save the solution
617 // internally in PHREEQC. The solution composition is saved after
618 // equilibration with the surface has been completed. The solution
619 // will not be saved for use in the next PHREEQC run.
620 os << "SAVE solution " << chemical_system_id + 1 << "\n";
621 }
622
623 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
624 if (!exchangers.empty())
625 {
626 os << "EXCHANGE " << chemical_system_id + 1 << "\n";
627 std::size_t const aqueous_solution_id =
628 dump->aqueous_solutions_prev.empty()
629 ? chemical_system_id + 1
630 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
631 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
632 for (auto const& exchanger : exchangers)
633 {
634 exchanger.print(os, chemical_system_id);
635 }
636 os << "SAVE solution " << chemical_system_id + 1 << "\n";
637 }
638
639 os << "END\n\n";
640 }
641
642 auto const& dump = phreeqc_io._dump;
643 if (dump)
644 {
645 dump->print(os, phreeqc_io._num_chemical_systems);
646 }
647
648 return os;
649}
650
652{
653 INFO("Phreeqc: Executing chemical calculation.");
654 if (RunFile(phreeqc_instance_id, _phreeqc_input_file.c_str()) != IPQ_OK)
655 {
656 OutputErrorString(phreeqc_instance_id);
657 OGS_FATAL(
658 "Failed in performing speciation calculation with the generated "
659 "phreeqc input file '{:s}'.",
661 }
662}
663
665{
666 auto const& basic_output_setups = _output->basic_output_setups;
667 auto const& phreeqc_result_file = basic_output_setups.output_file;
668 DBUG("Reading phreeqc results from file '{:s}'.", phreeqc_result_file);
669 std::ifstream in(phreeqc_result_file);
670
671 if (!in)
672 {
673 OGS_FATAL("Could not open phreeqc result file '{:s}'.",
674 phreeqc_result_file);
675 }
676
677 in >> *this;
678
679 if (!in)
680 {
681 OGS_FATAL("Error when reading phreeqc result file '{:s}'",
682 phreeqc_result_file);
683 }
684
685 in.close();
686}
687
688std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
689{
690 // Skip the headline
691 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
692
693 std::string line;
694 auto const& output = *phreeqc_io._output;
695 auto const& dropped_item_ids = output.dropped_item_ids;
696
697 auto const& surface = phreeqc_io._chemical_system->surface;
698 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
699
700 int const num_skipped_lines =
701 1 + (!surface.empty() ? 1 : 0) + (!exchangers.empty() ? 1 : 0);
702
703 auto& equilibrium_reactants =
704 phreeqc_io._chemical_system->equilibrium_reactants;
705 auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
706
707 for (std::size_t chemical_system_id = 0;
708 chemical_system_id < phreeqc_io._num_chemical_systems;
709 ++chemical_system_id)
710 {
711 // Skip equilibrium calculation result of initial solution
712 for (int i = 0; i < num_skipped_lines; ++i)
713 {
714 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
715 }
716
717 // Get calculation result of the solution after the reaction
718 if (!std::getline(in, line))
719 {
720 OGS_FATAL(
721 "Error when reading calculation result of Solution {:d} "
722 "after the reaction.",
723 chemical_system_id);
724 }
725
726 std::vector<double> accepted_items;
727 std::vector<std::string> items;
728 boost::trim_if(line, boost::is_any_of("\t "));
729 boost::algorithm::split(items, line, boost::is_any_of("\t "),
730 boost::token_compress_on);
731 for (int item_id = 0; item_id < static_cast<int>(items.size());
732 ++item_id)
733 {
734 if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
735 item_id) == dropped_item_ids.end())
736 {
737 double value;
738 try
739 {
740 value = std::stod(items[item_id]);
741 }
742 catch (const std::invalid_argument& e)
743 {
744 OGS_FATAL(
745 "Invalid argument. Could not convert string '{:s}' to "
746 "double for chemical system {:d}, column {:d}. "
747 "Exception '{:s}' was thrown.",
748 items[item_id], chemical_system_id + 1, item_id,
749 e.what());
750 }
751 catch (const std::out_of_range& e)
752 {
753 OGS_FATAL(
754 "Out of range error. Could not convert string "
755 "'{:s}' to double for chemical system {:d}, column "
756 "{:d}. Exception '{:s}' was thrown.",
757 items[item_id], chemical_system_id + 1, item_id,
758 e.what());
759 }
760 accepted_items.push_back(value);
761 }
762 }
763 assert(accepted_items.size() == output.accepted_items.size());
764
765 auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
766 auto& components = aqueous_solution->components;
767 auto& user_punch = phreeqc_io._user_punch;
768
769 GlobalIndexType const offset = aqueous_solution->pH->getRangeBegin();
770 GlobalIndexType const global_index = offset + chemical_system_id;
771
772 for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
773 ++item_id)
774 {
775 auto const& accepted_item = output.accepted_items[item_id];
776 auto const& item_name = accepted_item.name;
777
778 auto compare_by_name = [&item_name](auto const& item)
779 { return item.name == item_name; };
780
781 switch (accepted_item.item_type)
782 {
783 case ItemType::pH:
784 {
785 // Update pH value
787 *aqueous_solution->pH);
788 aqueous_solution->pH->set(
789 global_index, std::pow(10, -accepted_items[item_id]));
790 break;
791 }
792 case ItemType::pe:
793 {
794 // Update pe value
795 (*aqueous_solution->pe)[chemical_system_id] =
796 accepted_items[item_id];
797 break;
798 }
800 {
801 // Update component concentrations
802 auto const& component = BaseLib::findElementOrError(
803 components, compare_by_name,
804 [&]() {
805 OGS_FATAL("Could not find component '{:s}'.",
806 item_name);
807 });
809 *component.amount);
810 component.amount->set(global_index,
811 accepted_items[item_id]);
812 break;
813 }
815 {
816 // Update amounts of equilibrium reactant
817 auto const& equilibrium_reactant =
819 equilibrium_reactants, compare_by_name,
820 [&]()
821 {
822 OGS_FATAL(
823 "Could not find equilibrium reactant "
824 "'{:s}'",
825 item_name);
826 });
827 (*equilibrium_reactant.molality)[chemical_system_id] =
828 accepted_items[item_id];
829 break;
830 }
832 {
833 // Update amounts of kinetic reactants
834 auto const& kinetic_reactant = BaseLib::findElementOrError(
835 kinetic_reactants, compare_by_name,
836 [&]() {
837 OGS_FATAL("Could not find kinetic reactant '{:s}'.",
838 item_name);
839 });
840 (*kinetic_reactant.molality)[chemical_system_id] =
841 accepted_items[item_id];
842 break;
843 }
845 {
846 assert(user_punch);
847 auto const& secondary_variables =
848 user_punch->secondary_variables;
849 // Update values of secondary variables
850 auto const& secondary_variable =
852 secondary_variables, compare_by_name,
853 [&]() {
854 OGS_FATAL(
855 "Could not find secondary variable '{:s}'.",
856 item_name);
857 });
858 (*secondary_variable.value)[chemical_system_id] =
859 accepted_items[item_id];
860 break;
861 }
862 }
863 }
864 }
865
866 return in;
867}
868
869std::vector<std::string> const PhreeqcIO::getComponentList() const
870{
871 std::vector<std::string> component_names;
872 auto const& components = _chemical_system->aqueous_solution->components;
873 std::transform(components.begin(), components.end(),
874 std::back_inserter(component_names),
875 [](auto const& c) { return c.name; });
876
877 component_names.push_back("H");
878
879 return component_names;
880}
881
883 GlobalIndexType const& chemical_system_id,
884 MaterialPropertyLib::Medium const& medium,
885 ParameterLib::SpatialPosition const& pos, double const porosity,
886 double const t, double const dt)
887{
888 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
889 {
890 updateReactantVolumeFraction(kinetic_reactant, chemical_system_id,
891 medium, pos, porosity, t, dt);
892 }
893
894 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
895 {
896 updateReactantVolumeFraction(equilibrium_reactant, chemical_system_id,
897 medium, pos, porosity, t, dt);
898 }
899}
900
902 GlobalIndexType const& chemical_system_id,
903 MaterialPropertyLib::Medium const& medium,
904 double& porosity)
905{
906 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
907 {
908 setPorosityPostReaction(kinetic_reactant, chemical_system_id, medium,
909 porosity);
910 }
911
912 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
913 {
914 setPorosityPostReaction(equilibrium_reactant, chemical_system_id,
915 medium, porosity);
916 }
917}
918
920 std::size_t const ele_id,
921 std::vector<GlobalIndexType> const& chemical_system_indices)
922{
923 for (auto const& kinetic_reactant : _chemical_system->kinetic_reactants)
924 {
925 (*kinetic_reactant.mesh_prop_molality)[ele_id] =
926 averageReactantMolality(kinetic_reactant, chemical_system_indices);
927 }
928
929 for (auto const& equilibrium_reactant :
930 _chemical_system->equilibrium_reactants)
931 {
932 (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
933 averageReactantMolality(equilibrium_reactant,
934 chemical_system_indices);
935 }
936}
937} // namespace PhreeqcIOData
938} // namespace ChemistryLib
#define OGS_FATAL(...)
Definition Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the Mesh class.
std::vector< GlobalIndexType > chemical_system_index_map
double getConcentration(int const component_id, GlobalIndexType const chemical_system_id) const override
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)
void setAqueousSolutionsPrevFromDumpFile() override
std::unique_ptr< ChemicalSystem > _chemical_system
Definition PhreeqcIO.h:111
std::vector< ReactionRate > const _reaction_rates
Definition PhreeqcIO.h:112
std::unique_ptr< UserPunch > _user_punch
Definition PhreeqcIO.h:113
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
Definition PhreeqcIO.h:115
std::unique_ptr< Output > const _output
Definition PhreeqcIO.h:114
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::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
Definition Medium.cpp:33
Component const & component(std::size_t const &index) const
Definition Phase.cpp:33
ranges::range_reference_t< Range > findElementOrError(Range &range, std::predicate< ranges::range_reference_t< Range > > auto &&predicate, std::invocable auto error_callback)
Definition Algorithm.h:81
void initializeSiteMolality(Site &site, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Phase const &solid_phase, ParameterLib::SpatialPosition const &pos, double const t)
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)
Definition PhreeqcIO.cpp:86
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)
Definition PhreeqcIO.cpp:62
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::ostream & operator<<(std::ostream &os, PhreeqcIO const &phreeqc_io)
std::string specifyFileName(std::string const &project_file_name, std::string const &file_extension)
Definition Output.cpp:23
std::istream & operator>>(std::istream &in, PhreeqcIO &phreeqc_io)
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27