Add a migration physics to also handle the reaction term

refs idaholab#8
GiudGiud · Jan 13, 2025 · 11526f8 · 11526f8
1 parent dd70458
commit 11526f8
Show file tree

Hide file tree

Showing 2 changed files with 148 additions and 0 deletions.
diff --git a/include/physics/MultiSpeciesMigrationCG.h b/include/physics/MultiSpeciesMigrationCG.h
@@ -0,0 +1,27 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "MultiSpeciesDiffusionCG.h"
+
+/**
+ * Creates all the objects needed to solve diffusion equations for multiple species with a
+ * continuous Galerkin finite element discretization
+ */
+class MultiSpeciesMigrationCG : public MultiSpeciesDiffusionCG
+{
+public:
+  static InputParameters validParams();
+
+  MultiSpeciesMigrationCG(const InputParameters & parameters);
+
+private:
+  virtual void addFEKernels() override;
+};
diff --git a/src/physics/MultiSpeciesMigrationCG.C b/src/physics/MultiSpeciesMigrationCG.C
@@ -0,0 +1,121 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "MultiSpeciesMigrationCG.h"
+#include "MooseVariableBase.h"
+
+// Register the actions for the objects actually used
+registerMooseAction("MooseApp", MultiSpeciesMigrationCG, "add_kernel");
+registerMooseAction("MooseApp", MultiSpeciesMigrationCG, "add_bc");
+registerMooseAction("MooseApp", MultiSpeciesMigrationCG, "add_variable");
+registerMultiSpeciesDiffusionPhysicsBaseTasks("MooseApp", MultiSpeciesMigrationCG);
+
+InputParameters
+MultiSpeciesMigrationCG::validParams()
+{
+  InputParameters params = MultiSpeciesDiffusionCG::validParams();
+  params.addClassDescription(
+      "Discretizes a diffusion equation with the continuous Galerkin finite element method");
+
+  // Add reaction parameters
+  params.addParam<std::vector<std::vector<VariableName>>>(
+      "reacting_species", "For each species (outer indexing), the list of species they react with");
+  params.addParam<std::vector<std::vector<VariableName>>>(
+      "product_species", "For each species (outer indexing), for each reactant, the species being created");
+  params.addParam<std::vector<std::vector<MaterialPropertyName>>>(
+      "reaction_coefficients",
+      "For each species (outer indexing), the reaction coefficient for the reaction");
+
+  // Remove diffusion parameters for now: talk to PC
+  params.suppressParameter<std::vector<MooseFunctorName>>("diffusivity_functors");
+  params.suppressParameter<std::vector<MaterialPropertyName>>("diffusivity_matprops");
+
+  return params;
+}
+
+MultiSpeciesMigrationCG::MultiSpeciesMigrationCG(const InputParameters & parameters)
+  : MultiSpeciesDiffusionCG(parameters)
+{
+  checkTwoDVectorParamsSameLength<VariableName, MaterialPropertyName>("reacting_species",
+                                                                      "reaction_coefficients");
+  checkTwoDVectorParamsSameLength<VariableName, VariableName>("reacting_species",
+                                                              "product_species");
+}
+
+void
+MultiSpeciesMigrationCG::addFEKernels()
+{
+  MultiSpeciesDiffusionCG::addFEKernels();
+
+  if (!isParamValid("reacting_species"))
+    return;
+  const auto & reacting_species =
+      getParam<std::vector<std::vector<VariableName>>>("reacting_species");
+    const auto & product_species =
+        getParam<std::vector<std::vector<VariableName>>>("product_species");
+  const auto & reaction_coeffs =
+      getParam<std::vector<std::vector<MaterialPropertyName>>>("reaction_coefficients");
+
+  for (const auto s : index_range(_species_names))
+  {
+    const auto & var_name = _species_names[s];
+    // Reaction term
+    if (isParamValid("reacting_species"))
+    {
+      // Add a kernel for both directions
+      std::string kernel_type = "ADMatReactionFlexible";
+      InputParameters params = getFactory().getValidParams(kernel_type);
+      params.set<NonlinearVariableName>("variable") = var_name;
+      assignBlocks(params, _blocks);
+
+      // Double-indexed vectors arent initialized as expected
+      if (reacting_species.size() <= s)
+        continue;
+
+      for (const auto c : index_range(reacting_species[s]))
+      {
+        _console << s << " " << c << " " << reacting_species.size() << " " << reacting_species[s].size() << std::endl;
+        params.set<std::vector<VariableName>>("vs") = {var_name, reacting_species[s][c]};
+        params.set<MaterialPropertyName>("reaction_rate_name") = reaction_coeffs[s][c];
+        params.set<Real>("coeff") = -1;
+
+        getProblem().addKernel(kernel_type,
+                               prefix() + var_name + "_reaction_" + var_name + "_" +
+                                   reacting_species[s][c],
+                               params);
+
+        // only if the other reacting species is a nonlinear variable
+        if (nonlinearVariableExists(reacting_species[s][c], false))
+        {
+          params.set<NonlinearVariableName>("variable") = reacting_species[s][c];
+          params.set<std::vector<VariableName>>("vs") = {var_name, reacting_species[s][c]};
+
+          params.set<Real>("coeff") = -1;
+          getProblem().addKernel(kernel_type,
+                                 prefix() + var_name + "_reaction_" + reacting_species[s][c] + "_" +
+                                     var_name,
+                                 params);
+        }
+
+        // only if the target species is a nonlinear variable
+        if (nonlinearVariableExists(product_species[s][c], false))
+        {
+          params.set<NonlinearVariableName>("variable") = product_species[s][c];
+          params.set<std::vector<VariableName>>("vs") = {var_name, reacting_species[s][c]};
+
+          params.set<Real>("coeff") = 1;
+          getProblem().addKernel(kernel_type,
+                                 prefix() + var_name + "_production_" + var_name + "_" +
+                                     reacting_species[s][c],
+                                 params);
+        }
+      }
+    }
+  }
+}