feat: make output of H(0) and O(0) optional

ext/iphreeqc

@@ -1 +1 @@
-Subproject commit 38268b4aad03e6ce4755315f4cd690f007fa2720
+Subproject commit 2c8c311934cd54f0923d8b92dcccd9894825ce1b

src/Init/DiffusionInit.cpp

@@ -3,7 +3,6 @@
 // Rcpp include
 
 #include <map>
-#include <set>
 
 #include "DataStructures/Field.hpp"
 #include "InitialList.hpp"
@@ -15,7 +14,6 @@
 #include <Rinternals.h>
 #include <cstddef>
 #include <cstdint>
-#include <memory>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -26,9 +24,9 @@ namespace poet {
 enum SEXP_TYPE { SEXP_IS_LIST, SEXP_IS_VEC };
 
 const std::map<std::uint8_t, std::string> tug_side_mapping = {
-    {tug::BC_SIDE_RIGHT	, "E"},
-    {tug::BC_SIDE_LEFT	, "W"},
-    {tug::BC_SIDE_TOP	, "N"},
+    {tug::BC_SIDE_RIGHT, "E"},
+    {tug::BC_SIDE_LEFT, "W"},
+    {tug::BC_SIDE_TOP, "N"},
     {tug::BC_SIDE_BOTTOM, "S"}};
 
 static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
@@ -166,8 +164,8 @@ InitialList::resolveBoundaries(const Rcpp::List &boundaries_list,
       if (boundaries_list.containsElementNamed(side.second.c_str())) {
         const Rcpp::List mapping = boundaries_list[side.second];
 
-        const Rcpp::IntegerVector cells  = mapping["cell"];   // MDL 2024-06-13 
-        const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL 
+        const Rcpp::IntegerVector cells = mapping["cell"];    // MDL 2024-06-13
+        const Rcpp::IntegerVector values = mapping["sol_id"]; // MDL
         const Rcpp::CharacterVector type_str = mapping["type"];
 
         if (cells.size() != values.size()) {

src/Init/GridInit.cpp

@@ -6,10 +6,8 @@
 #include <Rcpp/Function.h>
 #include <Rcpp/vector/Matrix.h>
 #include <Rcpp/vector/instantiation.h>
-#include <cstdint>
 #include <fstream>
 #include <map>
-#include <memory>
 #include <regex>
 #include <sstream>
 #include <string>
@@ -123,7 +121,8 @@ static Rcpp::List expandGrid(const PhreeqcMatrix &pqc_mat,
   return Rcpp::Function("pqc_to_grid")(pqc_mat_R, grid_def);
 }
 
-PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
+PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input,
+                                       bool include_h0_o0) {
   // parse input values
   std::string script;
   std::string database;
@@ -180,7 +179,7 @@ PhreeqcMatrix InitialList::prepareGrid(const Rcpp::List &grid_input) {
     throw std::runtime_error("Grid size must be positive.");
   }
 
-  PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script);
+  PhreeqcMatrix pqc_mat = PhreeqcMatrix(database, script, include_h0_o0);
 
   this->transport_names = pqc_mat.getSolutionNames();
 

src/Init/InitialList.cpp

@@ -10,8 +10,9 @@
 #include <vector>
 
 namespace poet {
-void InitialList::initializeFromList(const Rcpp::List &setup) {
-  PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key]);
+void InitialList::initializeFromList(const Rcpp::List &setup,
+                                     bool include_h0_o0) {
+  PhreeqcMatrix phreeqc = prepareGrid(setup[grid_key], include_h0_o0);
   initDiffusion(setup[diffusion_key], phreeqc);
   initChemistry(setup[chemistry_key]);
 }
@@ -84,7 +85,8 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
   this->chem_hooks =
       Rcpp::as<Rcpp::List>(setup[static_cast<int>(ExportList::CHEM_HOOKS)]);
 
-  this->ai_surrogate_input_script = Rcpp::as<std::string>(setup[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)]);
+  this->ai_surrogate_input_script = Rcpp::as<std::string>(
+      setup[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)]);
 }
 
 Rcpp::List InitialList::exportList() {
@@ -132,7 +134,8 @@ Rcpp::List InitialList::exportList() {
   out[static_cast<int>(ExportList::CHEM_INTERP_SPECIES)] =
       Rcpp::wrap(this->interp_species);
   out[static_cast<int>(ExportList::CHEM_HOOKS)] = this->chem_hooks;
-  out[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)] = this->ai_surrogate_input_script;
+  out[static_cast<int>(ExportList::AI_SURROGATE_INPUT_SCRIPT)] =
+      this->ai_surrogate_input_script;
 
   return out;
 }

src/Init/InitialList.hpp

@@ -14,7 +14,6 @@
 #include <cstdint>
 
 #include <map>
-#include <memory>
 #include <string>
 #include <utility>
 #include <vector>
@@ -27,14 +26,14 @@ using TugType = double;
 
 class InitialList {
 public:
-  InitialList(RInside &R) : R(R){};
+  InitialList(RInside &R) : R(R) {};
 
-  void initializeFromList(const Rcpp::List &setup);
+  void initializeFromList(const Rcpp::List &setup, bool include_h0_o0 = false);
 
   void importList(const Rcpp::List &setup, bool minimal = false);
   Rcpp::List exportList();
 
-  Field getInitialGrid() const { return Field(this->initial_grid); }  
+  Field getInitialGrid() const { return Field(this->initial_grid); }
 
 private:
   RInside &R;
@@ -99,7 +98,7 @@ private:
 
   // std::unique_ptr<PhreeqcMatrix> pqc_mat;
 
-  PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input);
+  PhreeqcMatrix prepareGrid(const Rcpp::List &grid_input, bool include_h0_o0);
 
   std::uint8_t dim{0};
 
@@ -196,7 +195,7 @@ private:
   NamedVector<std::uint32_t> dht_species;
 
   NamedVector<std::uint32_t> interp_species;
-  
+
   // Path to R script that the user defines in the input file
   std::string ai_surrogate_input_script;
 

src/initializer.cpp

@@ -34,12 +34,15 @@ int main(int argc, char **argv) {
                "input script")
       ->default_val(false);
 
-  bool asRDS;
-  app.add_flag("-r, --rds", asRDS, "Save output as .rds")
-      ->default_val(false);
+  bool asRDS{false};
+  app.add_flag("-r, --rds", asRDS, "Save output as .rds")->default_val(false);
 
-  bool asQS;
-  app.add_flag("-q, --qs", asQS, "Save output as .qs")
+  bool asQS{false};
+  app.add_flag("-q, --qs", asQS, "Save output as .qs")->default_val(false);
+
+  bool includeH0O0;
+  app.add_flag("--include-h0-o0", includeH0O0,
+               "Include H(0) and O(0) in the output")
       ->default_val(false);
 
   CLI11_PARSE(app, argc, argv);
@@ -74,13 +77,13 @@ int main(int argc, char **argv) {
 
   // append the correct file extension
   if (asRDS) {
-      output_file += ".rds";
+    output_file += ".rds";
   } else if (asQS) {
-      output_file += ".qs";
+    output_file += ".qs";
   } else {
-      output_file += ".qs2";
+    output_file += ".qs2";
   }
-  
+
   // set working directory to the directory of the input script
   if (setwd) {
     const std::string dir_path = Rcpp::as<std::string>(