Merge pull request idaholab#17672 from YaqiWang/array_relaxation_17659

Fix array variable relaxation in Picard iteration

framework/include/loops/AllLocalDofIndicesThread.h

@@ -21,6 +21,7 @@ namespace libMesh
 class System;
 class DofMap;
 }
+class FEProblemBase;
 
 /**
  * Grab all the (possibly semi)local dof indices for the variables passed in, in the system passed
@@ -29,7 +30,7 @@ class DofMap;
 class AllLocalDofIndicesThread : public ParallelObject
 {
 public:
-  AllLocalDofIndicesThread(System & sys,
+  AllLocalDofIndicesThread(FEProblemBase & problem,
                            std::vector<std::string> vars,
                            bool include_semilocal = false);
   // Splitting Constructor
@@ -44,9 +45,8 @@ class AllLocalDofIndicesThread : public ParallelObject
   void dofIndicesSetUnion();
 
 protected:
-  System & _sys;
-  DofMap & _dof_map;
-  std::vector<std::string> _vars;
+  FEProblemBase & _problem;
+  System * _sys;
   std::vector<unsigned int> _var_numbers;
 
   /// Whether to include semilocal dof indices

framework/src/executioners/PicardSolve.C

@@ -165,8 +165,7 @@ PicardSolve::solve()
   if (_relax_factor != 1.0)
   {
     // Snag all of the local dof indices for all of these variables
-    System & libmesh_nl_system = _nl.system();
-    AllLocalDofIndicesThread aldit(libmesh_nl_system, _relaxed_vars);
+    AllLocalDofIndicesThread aldit(_problem, _relaxed_vars);
     ConstElemRange & elem_range = *_problem.mesh().getActiveLocalElementRange();
     Threads::parallel_reduce(elem_range, aldit);
 
@@ -178,8 +177,7 @@ PicardSolve::solve()
   if (_picard_self_relaxation_factor != 1.0)
   {
     // Snag all of the local dof indices for all of these variables
-    System & libmesh_nl_system = _nl.system();
-    AllLocalDofIndicesThread aldit(libmesh_nl_system, _picard_self_relaxed_variables);
+    AllLocalDofIndicesThread aldit(_problem, _picard_self_relaxed_variables);
     ConstElemRange & elem_range = *_problem.mesh().getActiveLocalElementRange();
     Threads::parallel_reduce(elem_range, aldit);
 

framework/src/executioners/Transient.C

@@ -20,7 +20,6 @@
 #include "Control.h"
 #include "TimePeriod.h"
 #include "MooseMesh.h"
-#include "AllLocalDofIndicesThread.h"
 #include "TimeIntegrator.h"
 #include "Console.h"
 

framework/src/kernels/ArrayDiffusion.C

@@ -17,9 +17,9 @@ InputParameters
 ArrayDiffusion::validParams()
 {
   InputParameters params = ArrayKernel::validParams();
-  params.addParam<MaterialPropertyName>("diffusion_coefficient",
-                                        "The name of the diffusivity, "
-                                        "can be scalar, vector, or matrix.");
+  params.addRequiredParam<MaterialPropertyName>("diffusion_coefficient",
+                                                "The name of the diffusivity, "
+                                                "can be scalar, vector, or matrix.");
   params.addClassDescription(
       "The array Laplacian operator ($-\\nabla \\cdot \\nabla u$), with the weak "
       "form of $(\\nabla \\phi_i, \\nabla u_h)$.");

framework/src/kernels/ArrayTimeDerivative.C

@@ -19,9 +19,9 @@ ArrayTimeDerivative::validParams()
   InputParameters params = ArrayTimeKernel::validParams();
   params.addClassDescription("Array time derivative operator with the weak form of $(\\psi_i, "
                              "\\frac{\\partial u_h}{\\partial t})$.");
-  params.addParam<MaterialPropertyName>("time_derivative_coefficient",
-                                        "The name of the time derivative coefficient. "
-                                        "Can be scalar, vector, or matrix");
+  params.addRequiredParam<MaterialPropertyName>("time_derivative_coefficient",
+                                                "The name of the time derivative coefficient. "
+                                                "Can be scalar, vector, or matrix");
   return params;
 }
 

framework/src/loops/AllLocalDofIndicesThread.C

@@ -11,6 +11,8 @@
 
 #include "FEProblem.h"
 #include "ParallelUniqueId.h"
+#include "NonlinearSystemBase.h"
+#include "MooseVariableFE.h"
 
 #include "libmesh/dof_map.h"
 #include "libmesh/threads.h"
@@ -21,27 +23,42 @@
 #include LIBMESH_INCLUDE_UNORDERED_SET
 LIBMESH_DEFINE_HASH_POINTERS
 
-AllLocalDofIndicesThread::AllLocalDofIndicesThread(System & sys,
+AllLocalDofIndicesThread::AllLocalDofIndicesThread(FEProblemBase & problem,
                                                    std::vector<std::string> vars,
                                                    bool include_semilocal)
-  : ParallelObject(sys.comm()),
-    _sys(sys),
-    _dof_map(sys.get_dof_map()),
-    _vars(vars),
+  : ParallelObject(problem.comm()),
+    _problem(problem),
+    _sys(nullptr),
     _include_semilocal(include_semilocal)
 {
-  _var_numbers.resize(_vars.size());
-  for (unsigned int i = 0; i < _vars.size(); i++)
-    _var_numbers[i] = _sys.variable_number(_vars[i]);
+  for (unsigned int i = 0; i < vars.size(); i++)
+  {
+    auto & var = _problem.getVariable(0, vars[i]);
+    if (_sys)
+    {
+      if (_sys != &var.sys().system())
+        mooseError("Variables passed in AllLocalDofIndicesThread must be all in the same system.");
+    }
+    else
+      _sys = &var.sys().system();
+
+    if (var.count() > 1) // array variable
+    {
+      const auto & array_var = _problem.getArrayVariable(0, vars[i]);
+      for (unsigned int p = 0; p < var.count(); ++p)
+        _var_numbers.push_back(_sys->variable_number(array_var.componentName(p)));
+    }
+    else
+      _var_numbers.push_back(_sys->variable_number(vars[i]));
+  }
 }
 
 // Splitting Constructor
 AllLocalDofIndicesThread::AllLocalDofIndicesThread(AllLocalDofIndicesThread & x,
                                                    Threads::split /*split*/)
-  : ParallelObject(x._sys.comm()),
+  : ParallelObject(x._problem.comm()),
+    _problem(x._problem),
     _sys(x._sys),
-    _dof_map(x._dof_map),
-    _vars(x._vars),
     _var_numbers(x._var_numbers),
     _include_semilocal(x._include_semilocal)
 {
@@ -53,17 +70,19 @@ AllLocalDofIndicesThread::operator()(const ConstElemRange & range)
   ParallelUniqueId puid;
   _tid = puid.id;
 
+  auto & dof_map = _sys->get_dof_map();
+
   for (const auto & elem : range)
   {
     std::vector<dof_id_type> dof_indices;
 
-    dof_id_type local_dof_begin = _dof_map.first_dof();
-    dof_id_type local_dof_end = _dof_map.end_dof();
+    dof_id_type local_dof_begin = dof_map.first_dof();
+    dof_id_type local_dof_end = dof_map.end_dof();
 
     // prepare variables
     for (unsigned int i = 0; i < _var_numbers.size(); i++)
     {
-      _dof_map.dof_indices(elem, dof_indices, _var_numbers[i]);
+      dof_map.dof_indices(elem, dof_indices, _var_numbers[i]);
       for (unsigned int j = 0; j < dof_indices.size(); j++)
       {
         dof_id_type dof = dof_indices[j];

framework/src/multiapps/TransientMultiApp.C

@@ -232,7 +232,7 @@ TransientMultiApp::solveStep(Real dt, Real target_time, bool auto_advance)
           transfer_old.close();
 
           // Snag all of the local dof indices for all of these variables
-          AllLocalDofIndicesThread aldit(libmesh_aux_system, _transferred_vars);
+          AllLocalDofIndicesThread aldit(problem, _transferred_vars);
           ConstElemRange & elem_range = *problem.mesh().getActiveLocalElementRange();
           Threads::parallel_reduce(elem_range, aldit);
 

framework/src/postprocessors/GreaterThanLessThanPostprocessor.C

@@ -66,7 +66,7 @@ GreaterThanLessThanPostprocessor::initialize()
 void
 GreaterThanLessThanPostprocessor::execute()
 {
-  AllLocalDofIndicesThread aldit(_fe_problem.getNonlinearSystemBase().system(), {_var.name()});
+  AllLocalDofIndicesThread aldit(_fe_problem, {_var.name()});
 
   if (_subdomain_restricted)
   {

framework/src/systems/NonlinearSystemBase.C

@@ -2680,7 +2680,7 @@ NonlinearSystemBase::computeJacobianInternal(const std::set<TagID> & tags)
 void
 NonlinearSystemBase::setVariableGlobalDoFs(const std::string & var_name)
 {
-  AllLocalDofIndicesThread aldit(_sys.system(), {var_name});
+  AllLocalDofIndicesThread aldit(_fe_problem, {var_name});
   ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
   Threads::parallel_reduce(elem_range, aldit);
 

framework/src/systems/SystemBase.C

@@ -223,25 +223,13 @@ SystemBase::getVariableBlocks(unsigned int var_number)
 void
 SystemBase::addVariableToZeroOnResidual(std::string var_name)
 {
-  unsigned int ncomp = getVariable(0, var_name).count();
-  if (ncomp > 1)
-    // need to push libMesh variable names for all components
-    for (unsigned int i = 0; i < ncomp; ++i)
-      _vars_to_be_zeroed_on_residual.push_back(_subproblem.arrayVariableComponent(var_name, i));
-  else
-    _vars_to_be_zeroed_on_residual.push_back(var_name);
+  _vars_to_be_zeroed_on_residual.push_back(var_name);
 }
 
 void
 SystemBase::addVariableToZeroOnJacobian(std::string var_name)
 {
-  unsigned int ncomp = getVariable(0, var_name).count();
-  if (ncomp > 1)
-    // need to push libMesh variable names for all components
-    for (unsigned int i = 0; i < ncomp; ++i)
-      _vars_to_be_zeroed_on_jacobian.push_back(_subproblem.arrayVariableComponent(var_name, i));
-  else
-    _vars_to_be_zeroed_on_jacobian.push_back(var_name);
+  _vars_to_be_zeroed_on_jacobian.push_back(var_name);
 }
 
 void
@@ -251,7 +239,11 @@ SystemBase::zeroVariables(std::vector<std::string> & vars_to_be_zeroed)
   {
     NumericVector<Number> & solution = this->solution();
 
-    AllLocalDofIndicesThread aldit(system(), vars_to_be_zeroed, true);
+    auto problem = dynamic_cast<FEProblemBase *>(&_subproblem);
+    if (!problem)
+      mooseError("System needs to be registered in FEProblemBase for using zeroVariables.");
+
+    AllLocalDofIndicesThread aldit(*problem, vars_to_be_zeroed, true);
     ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
     Threads::parallel_reduce(elem_range, aldit);
 

modules/contact/src/postprocessors/ContactDOFSetSize.C

@@ -54,7 +54,7 @@ ContactDOFSetSize::initialize()
 void
 ContactDOFSetSize::execute()
 {
-  AllLocalDofIndicesThread aldit(_fe_problem.getNonlinearSystemBase().system(), {_var.name()});
+  AllLocalDofIndicesThread aldit(_fe_problem, {_var.name()});
 
   // Get the element iterators corresponding to the subdomain id
   auto elem_begin = _mesh.active_local_subdomain_elements_begin(_subdomain_id);

test/include/auxkernels/ArrayQuotientAux.h

@@ -0,0 +1,26 @@
+//* 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 "AuxKernel.h"
+
+class ArrayQuotientAux : public ArrayAuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  ArrayQuotientAux(const InputParameters & parameters);
+
+protected:
+  virtual RealEigenVector computeValue() override;
+
+  const ArrayVariableValue & _numerator;
+  const ArrayVariableValue & _denominator;
+};

test/src/auxkernels/ArrayQuotientAux.C

@@ -0,0 +1,35 @@
+//* 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 "ArrayQuotientAux.h"
+
+registerMooseObject("MooseTestApp", ArrayQuotientAux);
+
+InputParameters
+ArrayQuotientAux::validParams()
+{
+  InputParameters params = ArrayAuxKernel::validParams();
+  params.addClassDescription("Divides two coupled variables.");
+  params.addCoupledVar("numerator", "The upstairs of the quotient variable");
+  params.addCoupledVar("denominator", "The downstairs of the quotient variable");
+  return params;
+}
+
+ArrayQuotientAux::ArrayQuotientAux(const InputParameters & parameters)
+  : ArrayAuxKernel(parameters),
+    _numerator(coupledArrayValue("numerator")),
+    _denominator(coupledArrayValue("denominator"))
+{
+}
+
+RealEigenVector
+ArrayQuotientAux::computeValue()
+{
+  return _numerator[_qp].cwiseQuotient(_denominator[_qp]);
+}