Merge pull request #241 from landinjm/fix_multiple_linear_solves

Fix multiple linear/nonlinear solves
prisms-center · Sep 16, 2024 · 6510515 · 6510515
2 parents e988b02 + 7f9b2cc
commit 6510515
Show file tree

Hide file tree

Showing 5 changed files with 220 additions and 368 deletions.
diff --git a/include/model_variables.h b/include/model_variables.h
@@ -34,7 +34,6 @@ class modelResidual
 struct variable_info
 {
   bool                                     is_scalar;
-  unsigned int                             variable_index;
   unsigned int                             global_var_index;
   dealii::EvaluationFlags::EvaluationFlags evaluation_flags;
   dealii::EvaluationFlags::EvaluationFlags residual_flags;

diff --git a/include/variableContainer.h b/include/variableContainer.h
@@ -3,11 +3,14 @@
 #ifndef VARIBLECONTAINER_H
 #define VARIBLECONTAINER_H
 
+#include <deal.II/base/exceptions.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 
+#include <boost/unordered_map.hpp>
+
 #include "userInputParameters.h"
 
 template <int dim, int degree, typename T>
@@ -19,15 +22,16 @@ class variableContainer
 
   // Standard contructor, used for most situations
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList,
-                    std::vector<variable_info>             _varChangeInfoList);
+                    const std::vector<variable_info>      &_varInfoList,
+                    const std::vector<variable_info>      &_varChangeInfoList);
+
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList);
+                    const std::vector<variable_info>      &_varInfoList);
   // Nonstandard constructor, used when only one index of "data" should be used,
   // use with care!
   variableContainer(const dealii::MatrixFree<dim, double> &data,
-                    std::vector<variable_info>             _varInfoList,
-                    unsigned int                           fixed_index);
+                    const std::vector<variable_info>      &_varInfoList,
+                    const unsigned int                    &fixed_index);
 
   // Methods to get the value/grad/hess in the residual method (this is how the
   // user gets these values in equations.h)
@@ -90,11 +94,6 @@ class variableContainer
   reinit_and_eval_change_in_solution(const vectorType &src,
                                      unsigned int      cell,
                                      unsigned int      var_being_solved);
-  void
-  reinit_and_eval_LHS(const vectorType               &src,
-                      const std::vector<vectorType *> solutionSet,
-                      unsigned int                    cell,
-                      unsigned int                    var_being_solved);
 
   // Only initialize the FEEvaluation object for each variable (used for
   // post-processing)
@@ -111,29 +110,34 @@ class variableContainer
   // The quadrature point index, a method to get the number of quadrature points
   // per cell, and a method to get the xyz coordinates for the quadrature point
   unsigned int q_point;
+
   unsigned int
-  get_num_q_points();
+  get_num_q_points() const;
+
   dealii::Point<dim, T>
-  get_q_point_location();
+  get_q_point_location() const;
 
 private:
-  // The number of variables
-  unsigned int num_var;
-
   // Vectors of the actual FEEvaluation objects for each active variable, split
   // into scalar variables and vector variables for type reasons
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>>   scalar_vars;
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>> vector_vars;
+  using scalar_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, 1, double>;
+  using vector_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, dim, double>;
+
+  boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>> scalar_vars_map;
+  boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>> vector_vars_map;
 
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>>
-    scalar_change_in_vars;
-  std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>>
-    vector_change_in_vars;
+  boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>>
+    scalar_change_in_vars_map;
+  boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>>
+    vector_change_in_vars_map;
 
   // Object containing some information about each variable (indices, whether
   // the val/grad/hess is needed, etc)
   std::vector<variable_info> varInfoList;
   std::vector<variable_info> varChangeInfoList;
+
+  // The number of variables
+  unsigned int num_var;
 };
 
 #endif
diff --git a/src/matrixfree/computeLHS.cc b/src/matrixfree/computeLHS.cc
@@ -57,7 +57,6 @@ MatrixFreePDE<dim, degree>::getLHS(
   for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
     {
       // Initialize, read DOFs, and set evaulation flags for each variable
-      // variable_list.reinit_and_eval_LHS(src,solutionSet,cell,currentFieldIndex);
       variable_list.reinit_and_eval(solutionSet, cell);
       variable_list.reinit_and_eval_change_in_solution(src, cell, currentFieldIndex);
 

diff --git a/src/userInputParameters/loadVariableAttributes.cc b/src/userInputParameters/loadVariableAttributes.cc
@@ -50,8 +50,6 @@ userInputParameters<dim>::loadVariableAttributes(
         }
     }
   varInfoListExplicitRHS.reserve(num_var_explicit_RHS);
-  unsigned int scalar_var_index = 0;
-  unsigned int vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -64,28 +62,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListExplicitRHS.push_back(varInfo);
     }
@@ -101,8 +80,6 @@ userInputParameters<dim>::loadVariableAttributes(
         }
     }
   varInfoListNonexplicitRHS.reserve(num_var_nonexplicit_RHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -115,28 +92,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListNonexplicitRHS.push_back(varInfo);
     }
@@ -153,8 +111,6 @@ userInputParameters<dim>::loadVariableAttributes(
     }
 
   varInfoListLHS.reserve(num_var_LHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -167,35 +123,14 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varInfoListLHS.push_back(varInfo);
     }
 
   varChangeInfoListLHS.reserve(num_var_LHS);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -209,37 +144,16 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       varChangeInfoListLHS.push_back(varInfo);
     }
 
   // Load variable information for postprocessing
   // First, the info list for the base field variables
   pp_baseVarInfoList.reserve(number_of_variables);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < number_of_variables; i++)
     {
       variable_info varInfo;
@@ -249,28 +163,9 @@ userInputParameters<dim>::loadVariableAttributes(
 
       varInfo.global_var_index = i;
 
-      !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
-        ? varInfo.var_needed = true
-        : varInfo.var_needed = false;
+      varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);
 
-      if (var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar = true;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = scalar_var_index;
-              scalar_var_index++;
-            }
-        }
-      else
-        {
-          varInfo.is_scalar = false;
-          if (varInfo.var_needed)
-            {
-              varInfo.variable_index = vector_var_index;
-              vector_var_index++;
-            }
-        }
+      varInfo.is_scalar = var_type[i] == SCALAR;
 
       pp_baseVarInfoList.push_back(varInfo);
     }
@@ -292,8 +187,6 @@ userInputParameters<dim>::loadVariableAttributes(
 
   // The info list for the postprocessing field variables
   pp_varInfoList.reserve(pp_number_of_variables);
-  scalar_var_index = 0;
-  vector_var_index = 0;
   for (unsigned int i = 0; i < pp_number_of_variables; i++)
     {
       variable_info varInfo;
@@ -303,18 +196,9 @@ userInputParameters<dim>::loadVariableAttributes(
         variable_attributes.equation_dependency_parser.eval_flags_residual_postprocess[i];
 
       varInfo.global_var_index = i;
-      if (pp_var_type[i] == SCALAR)
-        {
-          varInfo.is_scalar      = true;
-          varInfo.variable_index = scalar_var_index;
-          scalar_var_index++;
-        }
-      else
-        {
-          varInfo.is_scalar      = false;
-          varInfo.variable_index = vector_var_index;
-          vector_var_index++;
-        }
+
+      varInfo.is_scalar = pp_var_type[i] == SCALAR;
+
       pp_varInfoList.push_back(varInfo);
     }
 }