Merge pull request #25663 from parikshitbajpai/common-FE-FV_25661

FV variable support for ThermochimicaNodalData

framework/doc/content/source/interfaces/Coupleable.md

@@ -146,8 +146,7 @@ process. These functions can be found here:
 ## Writing directly to coupled variables
 
 Element- and nodal user objects as well AuxKernels may obtain a writable reference to a MOOSE field variable
-through the `Coupleable::writableVariable` function. The returned variable reference provides a `setNodalValue`
-method that can be used to set the nodal or elemental DOF value(s) of the variable.
+through the `Coupleable::writableVariable` function. The returned variable reference provides a `setDofValue` (for FE and FV variables) and `setNodalvalue` (only for FE variables) methods that can be used to set the nodal or elemental DOF value(s) of the variable.
 
 `Coupleable::writableVariable` enforces compatibility between the calling object type and the family of the
 requested variable. I.e. nodal user objects and AuxKernels may only obtain references to nodal variables, and
@@ -160,3 +159,6 @@ variable in the system may at most be written to by a single object on any given
 The user object and aux kernel thread loops check if an executed object has any writable variable references, and
 if so, will insert those variables into the aux solution vector. This obviates the need for using the
 [`ProjectionAux`](ProjectionAux.md) kernel.
+
+!alert warning
+`Coupleable::writableVariable` can let users write to both FE / FV from AuxKernels and UserObjects but one must exercise caution about whether Nodal or Elemental type AuxKernels / UOs are used as the quadrature would depend on this choice and might lead to segfault if a FV variable values are set using `setDofValue` function for non-zero values of `_qp` .

framework/include/interfaces/Coupleable.h

@@ -33,6 +33,10 @@ template <typename T>
 class DenseVector;
 }
 
+template <typename>
+class MooseVariableField;
+typedef MooseVariableField<Real> MooseWritableVariable;
+
 /**
  * Interface for objects that needs coupling capabilities
  *
@@ -492,16 +496,17 @@ class Coupleable
 
   /**
    * Returns a *writable* MooseVariable object for a nodal or elemental variable. Use
-   * var.setNodalValue(val[, idx]) in both cases (!) to set the solution DOF values. Only one object
-   * can obtain a writable reference in a simulation. Note that the written values will not ba
-   * available in the same system loop! E.g. values written using this API by a nodal AuxKernel will
-   * not be updated for other nodal AuxKernels during the same iteration over all nodes.
+   * var.setNodalValue(val[, idx]) in both cases (!) to set the solution DOF values. Only one
+   * object can obtain a writable reference in a simulation. Note that the written values will
+   * not ba available in the same system loop! E.g. values written using this API by a nodal
+   * AuxKernel will not be updated for other nodal AuxKernels during the same iteration over all
+   * nodes.
    * @param var_name Name of coupled variable
    * @param comp Component number for vector of coupled variables
-   * @return Reference to a MooseVariable for the coupled variable
+   * @return Reference to a MooseWritableVariable for the coupled variable
    * @see Kernel::value
    */
-  MooseVariable & writableVariable(const std::string & var_name, unsigned int comp = 0);
+  MooseWritableVariable & writableVariable(const std::string & var_name, unsigned int comp = 0);
 
   /**
    * Returns a *writable* reference to a coupled variable for writing to multiple
@@ -517,7 +522,7 @@ class Coupleable
   /**
    * Checks that the passed in variable is only accessed writable by one object in a given subdomain
    */
-  void checkWritableVar(MooseVariable * var);
+  void checkWritableVar(MooseWritableVariable * var);
 
   /**
    * Returns an old value from previous time step  of a coupled variable
@@ -1679,7 +1684,7 @@ class Coupleable
                                                    Moose::OLDER_SOLUTION_TAG};
 
   /// keep a set of allocated writable variable references to make sure only one object can obtain them per thread
-  std::vector<std::set<MooseVariable *>> _writable_coupled_variables;
+  std::vector<std::set<MooseWritableVariable *>> _writable_coupled_variables;
 };
 
 template <typename T>

framework/src/interfaces/Coupleable.C

@@ -847,10 +847,10 @@ Coupleable::coupledArrayValues(const std::string & var_name) const
   return coupledVectorHelper<const ArrayVariableValue *>(var_name, func);
 }
 
-MooseVariable &
+MooseWritableVariable &
 Coupleable::writableVariable(const std::string & var_name, unsigned int comp)
 {
-  auto * var = dynamic_cast<MooseVariable *>(getVar(var_name, comp));
+  auto * var = getVarHelper<MooseWritableVariable>(var_name, comp);
 
   const auto * aux = dynamic_cast<const AuxKernel *>(this);
   const auto * euo = dynamic_cast<const ElementUserObject *>(this);
@@ -920,7 +920,7 @@ Coupleable::writableCoupledValue(const std::string & var_name, unsigned int comp
 }
 
 void
-Coupleable::checkWritableVar(MooseVariable * var)
+Coupleable::checkWritableVar(MooseWritableVariable * var)
 {
   // check block restrictions for compatibility
   const auto * br = dynamic_cast<const BlockRestrictable *>(this);

modules/chemical_reactions/include/actions/ChemicalCompositionAction.h

@@ -68,6 +68,9 @@ class ChemicalCompositionAction : public Action
   /// Mass/amount unit
   std::string _munit;
 
+  /// Are the variables FV type
+  bool _is_fv;
+
   /// List of phases tracked by Thermochimica
   std::vector<std::string> _phases;
 

modules/chemical_reactions/include/userobjects/ThermochimicaNodalData.h

@@ -92,19 +92,19 @@ class ThermochimicaNodalData : public NodalUserObject
   ///@}
 
   /// Writable phase amount variables
-  std::vector<MooseVariable *> _ph;
+  std::vector<MooseWritableVariable *> _ph;
 
   /// Writable species amount variables
-  std::vector<MooseVariable *> _sp;
+  std::vector<MooseWritableVariable *> _sp;
 
   /// Writable vapour pressures for each element
-  std::vector<MooseVariable *> _vp;
+  std::vector<MooseWritableVariable *> _vp;
 
   /// Writable chemical potential variables for each element
-  std::vector<MooseVariable *> _el_pot;
+  std::vector<MooseWritableVariable *> _el_pot;
 
   /// Writable variable for molar amounts of each element in specified phase
-  std::vector<MooseVariable *> _el_ph;
+  std::vector<MooseWritableVariable *> _el_ph;
 
   /// Mass unit for output species
   const enum class OutputMassUnit { MOLES, FRACTION } _output_mass_unit;

modules/chemical_reactions/src/actions/ChemicalCompositionAction.C

@@ -13,7 +13,7 @@
 #include "MooseMesh.h"
 #include "MooseUtils.h"
 #include "MooseUtils.h"
-
+#include "AddVariableAction.h"
 #include "libmesh/string_to_enum.h"
 
 #ifdef THERMOCHIMICA_ENABLED
@@ -60,6 +60,7 @@ ChemicalCompositionAction::validParams()
   exec_enum = {EXEC_INITIAL, EXEC_TIMESTEP_END};
   params.addParam<ExecFlagEnum>(
       "execute_on", exec_enum, "When to execute the ThermochimicaNodalData UO");
+  params.addParam<bool>("is_fv", false, "Should the variables set up by action be of FV type");
 
   params.addParam<std::vector<std::string>>("output_phases", "List of phases to be output");
   params.addParam<std::vector<std::string>>(
@@ -87,6 +88,7 @@ ChemicalCompositionAction::ChemicalCompositionAction(const InputParameters & par
     _tunit(getParam<MooseEnum>("tunit")),
     _punit(getParam<MooseEnum>("punit")),
     _munit(getParam<MooseEnum>("munit")),
+    _is_fv(getParam<bool>("is_fv")),
     _phases(getParam<std::vector<std::string>>("output_phases")),
     _species(getParam<std::vector<std::string>>("output_species")),
     _output_mass_unit(getParam<MooseEnum>("output_species_unit")),
@@ -407,11 +409,14 @@ ChemicalCompositionAction::act()
   //
   if (_current_task == "add_aux_variable")
   {
-    const std::string aux_var_type = "MooseVariable";
+
+    auto aux_var_type = AddVariableAction::variableType(
+        FEType(Utility::string_to_enum<Order>(_problem->mesh().hasSecondOrderElements() ? "SECOND"
+                                                                                        : "FIRST"),
+               Utility::string_to_enum<FEFamily>("LAGRANGE")),
+        /* is_fv = */ _is_fv,
+        /* is_array = */ false);
     auto params = _factory.getValidParams(aux_var_type);
-    const bool second = _problem->mesh().hasSecondOrderElements();
-    params.set<MooseEnum>("order") = second ? "SECOND" : "FIRST";
-    params.set<MooseEnum>("family") = "LAGRANGE";
 
     for (const auto i : index_range(_elements))
       _problem->addAuxVariable(aux_var_type, _elements[i], params);

modules/chemical_reactions/src/userobjects/ThermochimicaNodalData.C

@@ -188,9 +188,9 @@ ThermochimicaNodalData::execute()
         mooseError("Failed to get index of phase '", _ph_names[i], "'");
       // Convert from 1-based (fortran) to 0-based (c++) indexing
       if (index - 1 < 0)
-        _ph[i]->setNodalValue(0.0, _qp);
+        _ph[i]->setDofValue(0.0, _qp);
       else
-        _ph[i]->setNodalValue(moles_phase[index - 1], _qp);
+        _ph[i]->setDofValue(moles_phase[index - 1], _qp);
     }
 
     auto db_phases = Thermochimica::getPhaseNamesSystem();
@@ -258,9 +258,9 @@ ThermochimicaNodalData::execute()
           _species_phase_pairs[i].second); // can we somehow use IDs instead of strings here?
 
       if (idbg == 0)
-        _sp[i]->setNodalValue(fraction, _qp);
+        _sp[i]->setDofValue(fraction, _qp);
       else if (idbg == 1)
-        _sp[i]->setNodalValue(0.0, _qp);
+        _sp[i]->setDofValue(0.0, _qp);
 #ifndef NDEBUG
       else
         mooseError("Failed to get phase speciation for phase '",
@@ -278,10 +278,10 @@ ThermochimicaNodalData::execute()
         auto [potential, idbg] = Thermochimica::getOutputChemPot(_element_potentials[i]);
 
         if (idbg == 0)
-          _el_pot[i]->setNodalValue(potential, _qp);
+          _el_pot[i]->setDofValue(potential, _qp);
         else if (idbg == 1)
           // element not present, just leave this at 0 for now
-          _el_pot[i]->setNodalValue(0.0, _qp);
+          _el_pot[i]->setDofValue(0.0, _qp);
         else if (idbg == -1)
           Moose::out << "getoutputchempot " << idbg << "\n";
       }
@@ -294,9 +294,9 @@ ThermochimicaNodalData::execute()
         libmesh_ignore(moles);
 
         if (idbg == 0)
-          _vp[i]->setNodalValue(fraction * pressure, _qp);
+          _vp[i]->setDofValue(fraction * pressure, _qp);
         else if (idbg == 1)
-          _vp[i]->setNodalValue(0.0, _qp);
+          _vp[i]->setDofValue(0.0, _qp);
 #ifndef NDEBUG
         else
           mooseError("Failed to get vapor pressure for phase '",
@@ -315,9 +315,9 @@ ThermochimicaNodalData::execute()
                                                                    _phase_element_pairs[i].first);
 
         if (idbg == 0)
-          _el_ph[i]->setNodalValue(moles, _qp);
+          _el_ph[i]->setDofValue(moles, _qp);
         else if (idbg == 1)
-          _el_ph[i]->setNodalValue(0.0, _qp);
+          _el_ph[i]->setDofValue(0.0, _qp);
 #ifndef NDEBUG
         else
           mooseError("Failed to get moles of element '",

modules/chemical_reactions/test/tests/thermochimica/MoRuPd.i

@@ -0,0 +1,56 @@
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 1
+    ny = 1
+  []
+[]
+
+[GlobalParams]
+  elements = 'Mo Ru Pd'
+  output_phases = 'BCCN HCPN'
+  output_species = 'HCPN:Pd'
+  output_element_potentials = 'mu:Pd'
+  output_vapor_pressures = 'vp:gas_ideal:Pd'
+  output_element_phases = 'ep:BCCN:Pd'
+[]
+
+[ChemicalComposition]
+  thermofile = Kaye_NobleMetals.dat
+  tunit = K
+  punit = atm
+  munit = moles
+  temperature = 2250
+  output_species_unit = mole_fraction
+[]
+
+[ICs]
+  [Mo]
+    type = FunctionIC
+    variable = Mo
+    function = '800*(1-x)+4.3*x'
+  []
+  [Ru]
+    type = FunctionIC
+    variable = Ru
+    function = '200*(1-x)+4.5*x'
+  []
+  [Pd]
+    type = ConstantIC
+    variable = Pd
+    value = 1.0e-8
+  []
+[]
+
+[Problem]
+  solve = false
+[]
+
+[Executioner]
+  type = Steady
+[]
+
+[Outputs]
+  exodus = true
+[]

modules/chemical_reactions/test/tests/thermochimica/tests

@@ -25,6 +25,18 @@
     required_submodule = 'contrib/thermochimica'
   []
 
+  [MoRuPd]
+    type = Exodiff
+    input = MoRuPd.i
+    exodiff = MoRuPd_out.e
+    design = 'ThermochimicaNodalData.md'
+    issues = '#25661'
+    requirement = 'The system shall be able to use graciously handle a missing element in Thermochimica'
+    max_threads = 1
+    rel_err = 1e-3 # this is the same error thermochimica uses for its internal tests
+    required_submodule = 'contrib/thermochimica'
+  []
+
   [csv_ic]
     type = Exodiff
     input = csv_ic.i

test/include/auxkernels/MultipleUpdateAux.h

@@ -31,8 +31,8 @@ class MultipleUpdateAux : public AuxKernel
   const bool _deprecated;
 
   /// current API
-  MooseVariable * _var1;
-  MooseVariable * _var2;
+  MooseWritableVariable * _var1;
+  MooseWritableVariable * _var2;
 
   /// deprectated API
   VariableValue * _dvar1;

test/include/auxkernels/MultipleUpdateElemAux.h

@@ -28,7 +28,7 @@ class MultipleUpdateElemAux : public AuxKernel
   virtual void computeVarValues(std::vector<Real> & values);
 
   unsigned int _n_vars;
-  std::vector<MooseVariable *> _vars;
+  std::vector<MooseWritableVariable *> _vars;
 
   const bool _use_compute_value;
 };

test/include/kernels/MultipleUpdateErrorKernel.h

@@ -30,8 +30,8 @@ class MultipleUpdateErrorKernel : public Kernel
   const bool _deprecated;
 
   /// current API
-  MooseVariable * _var1;
-  MooseVariable * _var2;
+  MooseWritableVariable * _var1;
+  MooseWritableVariable * _var2;
 
   /// deprectated API
   VariableValue * _dvar1;

test/include/userobjects/MultiUpdateElementalUO.h

@@ -27,5 +27,5 @@ class MultiUpdateElementalUO : public ElementUserObject
   virtual void threadJoin(const UserObject &) override {}
 
 protected:
-  MooseVariable & _v;
+  MooseWritableVariable & _v;
 };

test/include/userobjects/MultiUpdateNodalUO.h

@@ -27,5 +27,5 @@ class MultiUpdateNodalUO : public NodalUserObject
   virtual void threadJoin(const UserObject &) override {}
 
 protected:
-  MooseVariable & _v;
+  MooseWritableVariable & _v;
 };

test/src/auxkernels/MultipleUpdateAux.C

@@ -52,8 +52,14 @@ MultipleUpdateAux::computeValue()
   }
   else
   {
-    _var1->setNodalValue(_nl_u[_qp] + 10.0, _qp);
-    _var2->setNodalValue(_nl_u[_qp] + 200.0, _qp);
+    /*
+      For NodalKernels the index _qp is always 0 and the computeValue method is executed on each
+      node but when using ElementalKernels the computeValue method is executed on each quadrature
+      point of an element. For this reason, in multi_update_fv_test.i input file, the quadrature is
+      set to Constant order since for FV variables there's only one DOF value locally.
+    */
+    _var1->setDofValue(_nl_u[_qp] + 10.0, _qp);
+    _var2->setDofValue(_nl_u[_qp] + 200.0, _qp);
   }
   return -3.33;
 }