Thermochimica API enhancements

modules/chemical_reactions/include/actions/ChemicalCompositionAction.h

@@ -49,4 +49,7 @@ class ChemicalCompositionAction : public Action
 
   /// List of element chemical potentials to be extracted from Thermochimica
   std::vector<std::string> _element_potentials;
+
+  /// List of gas phase species to extract vapor pressures from Thermochimica
+  std::vector<std::string> _vapor_pressures;
 };

modules/chemical_reactions/include/auxkernels/ThermochimicaAux.h

@@ -29,23 +29,22 @@ class ThermochimicaAux : public AuxKernel
   /// Total number of phases to output
   const std::size_t _n_phases;
   /// Writable phase amount variables
-  std::vector<VariableValue *> _ph;
-  /// Phase names
-  std::vector<std::string> _ph_name;
+  std::vector<MooseVariable *> _ph;
 
   /// Total number of species to output
   const std::size_t _n_species;
   /// Writable species amount variables
-  std::vector<VariableValue *> _sp;
-  /// Corresponding phase names
-  std::vector<std::string> _sp_phase_name;
-  /// Corresponding species names
-  std::vector<std::string> _sp_species_name;
+  std::vector<MooseVariable *> _sp;
+
+  /// Total number of vapor species
+  const std::size_t _n_vapor_species;
+  /// Writable vapour pressures for each element
+  std::vector<MooseVariable *> _vapor_pressures;
 
   /// Total number of elements to output
   const std::size_t _n_elements;
   /// Writable chemical potential variables for each element
-  std::vector<VariableValue *> _el_pot;
+  std::vector<MooseVariable *> _el_pot;
 
 private:
 #ifdef THERMOCHIMICA_ENABLED

modules/chemical_reactions/include/userobjects/ThermochimicaNodalData.h

@@ -51,14 +51,15 @@ class ThermochimicaNodalData : public NodalUserObject
     std::vector<double> _chemical_potential;
     std::vector<double> _mol_fraction;
     std::vector<Real> _species_fractions;
+    std::vector<Real> _vapor_pressures;
     std::vector<int> _phase_indices;
     std::vector<double> _element_potential_for_output;
   };
 
   const Data & getNodalData(dof_id_type node_id) const;
 
 protected:
-  Real _pressure;
+  const VariableValue & _pressure;
   const VariableValue & _temperature;
 
   // re-initialization data
@@ -67,19 +68,28 @@ class ThermochimicaNodalData : public NodalUserObject
   const std::size_t _n_phases;
   const std::size_t _n_species;
   const std::size_t _n_elements;
+  const std::size_t _n_vapor_species;
 
   std::vector<const VariableValue *> _el;
   std::vector<std::string> _el_name;
+  std::vector<unsigned int> _el_id;
+
+  std::pair<int, int> _db_num_phases;
+  std::vector<std::string> _db_phase_names;
+  std::vector<std::vector<std::string>> _db_species_names;
 
   std::vector<std::string> _ph_name;
   std::vector<std::string> _sp_phase_name;
   std::vector<std::string> _sp_species_name;
+  std::vector<std::string> _vapor_phase_name;
+  std::vector<std::string> _vapor_species_name;
 
   /// Nodal data (TODO: investigate writing directly to AuxVariables)
   std::unordered_map<dof_id_type, Data> _data;
 
   ///@{ Element chemical potential output
   const bool _output_element_potential;
+  const bool _output_vapor_pressures;
   std::vector<std::string> _element_potentials;
   ///@}
 };

modules/chemical_reactions/src/actions/ChemicalCompositionAction.C

@@ -40,21 +40,23 @@ ChemicalCompositionAction::validParams()
   params.addParam<FileName>("initial_values", "The CSV file name with initial conditions.");
   params.addParam<FileName>("thermofile", "Thermodynamics model file");
 
-  MooseEnum tUnit("K C F R", "K");
-  params.addParam<MooseEnum>("tunit", tUnit, "Temperature Unit");
-  MooseEnum pUnit("atm psi bar Pa kPa", "atm");
-  params.addParam<MooseEnum>("punit", pUnit, "Pressure Unit");
+  MooseEnum tUnit("K C F R");
+  params.addRequiredParam<MooseEnum>("tunit", tUnit, "Temperature Unit");
+  MooseEnum pUnit("atm psi bar Pa kPa");
+  params.addRequiredParam<MooseEnum>("punit", pUnit, "Pressure Unit");
   MooseEnum mUnit(
-      "mole_fraction atom_fraction atoms moles gram-atoms mass_fraction kilograms grams pounds",
-      "moles");
-  params.addParam<MooseEnum>("munit", mUnit, "Mass Unit");
+      "mole_fraction atom_fraction atoms moles gram-atoms mass_fraction kilograms grams pounds");
+  params.addRequiredParam<MooseEnum>("munit", mUnit, "Mass Unit");
 
   params.addParam<std::vector<std::string>>("output_phases", "List of phases to be output");
   params.addParam<std::vector<std::string>>(
       "output_species", "List species for which concentration in the phases is needed");
   params.addParam<std::vector<std::string>>(
       "element_potentials",
       "List of chemical elements for which chemical potentials are requested");
+  params.addParam<std::vector<std::string>>(
+      "output_vapor_pressures",
+      "List of gas phase species for which vapor pressures are requested");
   return params;
 }
 
@@ -66,7 +68,8 @@ ChemicalCompositionAction::ChemicalCompositionAction(const InputParameters & par
     _munit(getParam<MooseEnum>("munit")),
     _phases(getParam<std::vector<std::string>>("output_phases")),
     _species(getParam<std::vector<std::string>>("output_species")),
-    _element_potentials(getParam<std::vector<std::string>>("element_potentials"))
+    _element_potentials(getParam<std::vector<std::string>>("element_potentials")),
+    _vapor_pressures(getParam<std::vector<std::string>>("output_vapor_pressures"))
 {
   ThermochimicaUtils::checkLibraryAvailability(*this);
 
@@ -99,6 +102,9 @@ ChemicalCompositionAction::act()
 
     for (const auto i : index_range(_element_potentials))
       _problem->addAuxVariable(aux_var_type, _element_potentials[i], params);
+
+    for (const auto i : index_range(_vapor_pressures))
+      _problem->addAuxVariable(aux_var_type, _vapor_pressures[i], params);
   }
 
   //
@@ -180,6 +186,13 @@ ChemicalCompositionAction::act()
       params.set<AuxVariableName>("variable") = ker_name;
       _problem->addAuxKernel("SelfAux", ker_name, params);
     }
+
+    for (const auto i : index_range(_vapor_pressures))
+    {
+      const std::string ker_name = _vapor_pressures[i];
+      params.set<AuxVariableName>("variable") = ker_name;
+      _problem->addAuxKernel("SelfAux", ker_name, params);
+    }
   }
 
 #endif

modules/chemical_reactions/src/auxkernels/ThermochimicaAux.C

@@ -31,6 +31,7 @@ ThermochimicaAux::validParams()
   params.addCoupledVar("output_species", "Amounts of species to be output");
   params.addRequiredParam<UserObjectName>("thermo_nodal_data_uo", "Name of the user object");
   params.addCoupledVar("element_potentials", "Chemical potentials of elements");
+  params.addCoupledVar("output_vapor_pressures", "Vapour pressures of species to be output");
 
   ThermochimicaUtils::addClassDescription(
       params,
@@ -44,11 +45,10 @@ ThermochimicaAux::ThermochimicaAux(const InputParameters & parameters)
   : AuxKernel(parameters),
     _n_phases(coupledComponents("output_phases")),
     _ph(_n_phases),
-    _ph_name(_n_phases),
     _n_species(coupledComponents("output_species")),
     _sp(_n_species),
-    _sp_phase_name(_n_species),
-    _sp_species_name(_n_species),
+    _n_vapor_species(coupledComponents("output_vapor_pressures")),
+    _vapor_pressures(_n_vapor_species),
     _n_elements(coupledComponents("element_potentials")),
     _el_pot(_n_elements)
 #ifdef THERMOCHIMICA_ENABLED
@@ -65,22 +65,16 @@ ThermochimicaAux::ThermochimicaAux(const InputParameters & parameters)
       mooseError("All variables coupled in ThermochimicaAux must be of first order Lagrange type.");
 
   for (const auto i : make_range(_n_phases))
-  {
-    _ph[i] = &writableCoupledValue("output_phases", i);
-    _ph_name[i] = getVar("output_phases", i)->name();
-  }
+    _ph[i] = &writableVariable("output_phases", i);
 
   for (const auto i : make_range(_n_species))
-  {
-    _sp[i] = &writableCoupledValue("output_species", i);
-    auto species_var_name = getVar("output_species", i)->name();
-    int semicolon = species_var_name.find(";");
-    _sp_phase_name[i] = species_var_name.substr(0, semicolon);
-    _sp_species_name[i] = species_var_name.substr(semicolon + 1);
-  }
+    _sp[i] = &writableVariable("output_species", i);
+
+  for (const auto i : make_range(_n_vapor_species))
+    _vapor_pressures[i] = &writableVariable("output_vapor_pressures", i);
 
   for (const auto i : make_range(_n_elements))
-    _el_pot[i] = &writableCoupledValue("element_potentials", i);
+    _el_pot[i] = &writableVariable("element_potentials", i);
 }
 
 Real
@@ -95,21 +89,28 @@ ThermochimicaAux::computeValue()
   // Save requested phase data into coupled aux variables
   for (const auto i : make_range(_n_phases))
     if (data._phase_indices[i] < 0)
-      (*_ph[i])[_qp] = 0;
+      _ph[i]->setNodalValue(0, _qp);
     else
     {
-      (*_ph[i])[_qp] = data._moles_phase[data._phase_indices[i]];
+      _ph[i]->setNodalValue(data._moles_phase[data._phase_indices[i]], _qp);
       n_active_phases += 1.0;
     }
 
   // Save requested species data into coupled aux variables
   for (unsigned int i = 0; i < _n_species; i++)
-    (*_sp[i])[_qp] = data._species_fractions[i];
+    _sp[i]->setNodalValue(data._species_fractions[i], _qp);
+
+  // Save requested vapor pressures into coupled aux variables
+  mooseAssert(_vapor_pressures.size() == data._vapor_pressures.size(),
+              "Output vapor pressures: Inconsistent sizes.");
+  for (const auto i : make_range(_n_vapor_species))
+    _vapor_pressures[i]->setNodalValue(data._vapor_pressures[i], _qp);
 
   // Save requested element potentials into coupled aux variables
-  mooseAssert(_el_pot.size() == data._element_potential.size(), "Inconsistent sizes.");
+  mooseAssert(_el_pot.size() == data._element_potential_for_output.size(),
+              "Output element potentials: Inconsistent sizes.");
   for (const auto i : make_range(_n_elements))
-    (*_el_pot[i])[_qp] = data._element_potential[i];
+    _el_pot[i]->setNodalValue(data._element_potential_for_output[i], _qp);
 #endif
 
   return n_active_phases;