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. Normally, this happens in
445 // the first time step.
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>::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 if (!aqueous_solutions_prev.empty())
525 {
526 os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
527 }
528 }
529
530 os << "USE solution none\n";
531 os << "END\n\n";
532
533 os << "USE solution " << chemical_system_id + 1 << "\n\n";
534
535 auto const& equilibrium_reactants =
536 phreeqc_io._chemical_system->equilibrium_reactants;
537 if (!equilibrium_reactants.empty() || fixing_pe)
538 {
539 os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
540 for (auto const& equilibrium_reactant : equilibrium_reactants)
541 {
542 equilibrium_reactant.print(os, chemical_system_id);
543 }
544 fixing_pe
545 ? os << "Fix_pe "
546 << -phreeqc_io._chemical_system->aqueous_solution->pe0
547 << " O2(g)\n\n"
548 : os << "\n";
549 }
550
551 auto const& kinetic_reactants =
552 phreeqc_io._chemical_system->kinetic_reactants;
553 if (!kinetic_reactants.empty())
554 {
555 os << "KINETICS " << chemical_system_id + 1 << "\n";
556 for (auto const& kinetic_reactant : kinetic_reactants)
557 {
558 kinetic_reactant.print(os, chemical_system_id);
559 }
560 os << "-steps " << phreeqc_io._dt << "\n\n";
561 }
562
563 auto const& surface = phreeqc_io._chemical_system->surface;
564 if (!surface.empty())
565 {
566 os << "SURFACE " << chemical_system_id + 1 << "\n";
567 std::size_t aqueous_solution_id =
568 dump->aqueous_solutions_prev.empty()
569 ? chemical_system_id + 1
570 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
571 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
572
573 // print unit
574 if (std::holds_alternative<DensityBasedSurfaceSite>(
575 surface.front()))
576 {
577 os << "-sites_units density\n";
578 }
579 else
580 {
581 os << "-sites_units absolute\n";
582 }
583
584 for (auto const& surface_site : surface)
585 {
586 std::visit(
587 overloaded{[&os](DensityBasedSurfaceSite const& s)
588 {
589 os << s.name << " " << s.site_density << " "
590 << s.specific_surface_area << " "
591 << s.mass << "\n";
592 },
593 [&os, chemical_system_id](
594 MoleBasedSurfaceSite const& s) {
595 os << s.name << " "
596 << (*s.molality)[chemical_system_id]
597 << "\n";
598 }},
599 surface_site);
600 }
601
602 // overlook the effect of the buildup of charges onto the surface
603 if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
604 {
605 os << "-no_edl\n";
606 }
607
608 os << "SAVE solution " << chemical_system_id + 1 << "\n";
609 }
610
611 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
612 if (!exchangers.empty())
613 {
614 os << "EXCHANGE " << chemical_system_id + 1 << "\n";
615 std::size_t const aqueous_solution_id =
616 dump->aqueous_solutions_prev.empty()
617 ? chemical_system_id + 1
618 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
619 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
620 for (auto const& exchanger : exchangers)
621 {
622 exchanger.print(os, chemical_system_id);
623 }
624 os << "SAVE solution " << chemical_system_id + 1 << "\n";
625 }
626
627 os << "END\n\n";
628 }
629
630 auto const& dump = phreeqc_io._dump;
631 if (dump)
632 {
633 dump->print(os, phreeqc_io._num_chemical_systems);
634 }
635
636 return os;
637}
638
640{
641 INFO("Phreeqc: Executing chemical calculation.");
642 if (RunFile(phreeqc_instance_id, _phreeqc_input_file.c_str()) != IPQ_OK)
643 {
644 OutputErrorString(phreeqc_instance_id);
645 OGS_FATAL(
646 "Failed in performing speciation calculation with the generated "
647 "phreeqc input file '{:s}'.",
649 }
650}
651
653{
654 auto const& basic_output_setups = _output->basic_output_setups;
655 auto const& phreeqc_result_file = basic_output_setups.output_file;
656 DBUG("Reading phreeqc results from file '{:s}'.", phreeqc_result_file);
657 std::ifstream in(phreeqc_result_file);
658
659 if (!in)
660 {
661 OGS_FATAL("Could not open phreeqc result file '{:s}'.",
662 phreeqc_result_file);
663 }
664
665 in >> *this;
666
667 if (!in)
668 {
669 OGS_FATAL("Error when reading phreeqc result file '{:s}'",
670 phreeqc_result_file);
671 }
672
673 in.close();
674}
675
676std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
677{
678 // Skip the headline
679 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
680
681 std::string line;
682 auto const& output = *phreeqc_io._output;
683 auto const& dropped_item_ids = output.dropped_item_ids;
684
685 auto const& surface = phreeqc_io._chemical_system->surface;
686 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
687
688 int const num_skipped_lines =
689 1 + (!surface.empty() ? 1 : 0) + (!exchangers.empty() ? 1 : 0);
690
691 auto& equilibrium_reactants =
692 phreeqc_io._chemical_system->equilibrium_reactants;
693 auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
694
695 for (std::size_t chemical_system_id = 0;
696 chemical_system_id < phreeqc_io._num_chemical_systems;
697 ++chemical_system_id)
698 {
699 // Skip equilibrium calculation result of initial solution
700 for (int i = 0; i < num_skipped_lines; ++i)
701 {
702 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
703 }
704
705 // Get calculation result of the solution after the reaction
706 if (!std::getline(in, line))
707 {
708 OGS_FATAL(
709 "Error when reading calculation result of Solution {:d} "
710 "after the reaction.",
711 chemical_system_id);
712 }
713
714 std::vector<double> accepted_items;
715 std::vector<std::string> items;
716 boost::trim_if(line, boost::is_any_of("\t "));
717 boost::algorithm::split(items, line, boost::is_any_of("\t "),
718 boost::token_compress_on);
719 for (int item_id = 0; item_id < static_cast<int>(items.size());
720 ++item_id)
721 {
722 if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
723 item_id) == dropped_item_ids.end())
724 {
725 double value;
726 try
727 {
728 value = std::stod(items[item_id]);
729 }
730 catch (const std::invalid_argument& e)
731 {
732 OGS_FATAL(
733 "Invalid argument. Could not convert string '{:s}' to "
734 "double for chemical system {:d}, column {:d}. "
735 "Exception '{:s}' was thrown.",
736 items[item_id], chemical_system_id + 1, item_id,
737 e.what());
738 }
739 catch (const std::out_of_range& e)
740 {
741 OGS_FATAL(
742 "Out of range error. Could not convert string "
743 "'{:s}' to double for chemical system {:d}, column "
744 "{:d}. Exception '{:s}' was thrown.",
745 items[item_id], chemical_system_id + 1, item_id,
746 e.what());
747 }
748 accepted_items.push_back(value);
749 }
750 }
751 assert(accepted_items.size() == output.accepted_items.size());
752
753 auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
754 auto& components = aqueous_solution->components;
755 auto& user_punch = phreeqc_io._user_punch;
756
757 GlobalIndexType const offset = aqueous_solution->pH->getRangeBegin();
758 GlobalIndexType const global_index = offset + chemical_system_id;
759
760 for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
761 ++item_id)
762 {
763 auto const& accepted_item = output.accepted_items[item_id];
764 auto const& item_name = accepted_item.name;
765
766 auto compare_by_name = [&item_name](auto const& item)
767 { return item.name == item_name; };
768
769 switch (accepted_item.item_type)
770 {
771 case ItemType::pH:
772 {
773 // Update pH value
775 *aqueous_solution->pH);
776 aqueous_solution->pH->set(
777 global_index, std::pow(10, -accepted_items[item_id]));
778 break;
779 }
780 case ItemType::pe:
781 {
782 // Update pe value
783 (*aqueous_solution->pe)[chemical_system_id] =
784 accepted_items[item_id];
785 break;
786 }
788 {
789 // Update component concentrations
790 auto const& component = BaseLib::findElementOrError(
791 components, compare_by_name,
792 [&]() {
793 OGS_FATAL("Could not find component '{:s}'.",
794 item_name);
795 });
797 *component.amount);
798 component.amount->set(global_index,
799 accepted_items[item_id]);
800 break;
801 }
803 {
804 // Update amounts of equilibrium reactant
805 auto const& equilibrium_reactant =
807 equilibrium_reactants, compare_by_name,
808 [&]()
809 {
810 OGS_FATAL(
811 "Could not find equilibrium reactant "
812 "'{:s}'",
813 item_name);
814 });
815 (*equilibrium_reactant.molality)[chemical_system_id] =
816 accepted_items[item_id];
817 break;
818 }
820 {
821 // Update amounts of kinetic reactants
822 auto const& kinetic_reactant = BaseLib::findElementOrError(
823 kinetic_reactants, compare_by_name,
824 [&]() {
825 OGS_FATAL("Could not find kinetic reactant '{:s}'.",
826 item_name);
827 });
828 (*kinetic_reactant.molality)[chemical_system_id] =
829 accepted_items[item_id];
830 break;
831 }
833 {
834 assert(user_punch);
835 auto const& secondary_variables =
836 user_punch->secondary_variables;
837 // Update values of secondary variables
838 auto const& secondary_variable =
840 secondary_variables, compare_by_name,
841 [&]() {
842 OGS_FATAL(
843 "Could not find secondary variable '{:s}'.",
844 item_name);
845 });
846 (*secondary_variable.value)[chemical_system_id] =
847 accepted_items[item_id];
848 break;
849 }
850 }
851 }
852 }
853
854 return in;
855}
856
857std::vector<std::string> const PhreeqcIO::getComponentList() const
858{
859 std::vector<std::string> component_names;
860 auto const& components = _chemical_system->aqueous_solution->components;
861 std::transform(components.begin(), components.end(),
862 std::back_inserter(component_names),
863 [](auto const& c) { return c.name; });
864
865 component_names.push_back("H");
866
867 return component_names;
868}
869
871 GlobalIndexType const& chemical_system_id,
872 MaterialPropertyLib::Medium const& medium,
873 ParameterLib::SpatialPosition const& pos, double const porosity,
874 double const t, double const dt)
875{
876 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
877 {
878 updateReactantVolumeFraction(kinetic_reactant, chemical_system_id,
879 medium, pos, porosity, t, dt);
880 }
881
882 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
883 {
884 updateReactantVolumeFraction(equilibrium_reactant, chemical_system_id,
885 medium, pos, porosity, t, dt);
886 }
887}
888
890 GlobalIndexType const& chemical_system_id,
891 MaterialPropertyLib::Medium const& medium,
892 double& porosity)
893{
894 for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
895 {
896 setPorosityPostReaction(kinetic_reactant, chemical_system_id, medium,
897 porosity);
898 }
899
900 for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
901 {
902 setPorosityPostReaction(equilibrium_reactant, chemical_system_id,
903 medium, porosity);
904 }
905}
906
908 std::size_t const ele_id,
909 std::vector<GlobalIndexType> const& chemical_system_indices)
910{
911 for (auto const& kinetic_reactant : _chemical_system->kinetic_reactants)
912 {
913 (*kinetic_reactant.mesh_prop_molality)[ele_id] =
914 averageReactantMolality(kinetic_reactant, chemical_system_indices);
915 }
916
917 for (auto const& equilibrium_reactant :
918 _chemical_system->equilibrium_reactants)
919 {
920 (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
921 averageReactantMolality(equilibrium_reactant,
922 chemical_system_indices);
923 }
924}
925} // namespace PhreeqcIOData
926} // 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:76
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