Create utility to generate map from global dof index to derivative

Also move ADOffset to be ADUtils

Closes #15450

framework/include/materials/MaterialProperty.h

@@ -109,6 +109,8 @@ template <typename T, bool is_ad>
 class MaterialPropertyBase : public PropertyValue
 {
 public:
+  typedef MooseADWrapper<T, is_ad> value_type;
+
   /// Explicitly declare a public constructor because we made the copy constructor private
   MaterialPropertyBase() : PropertyValue()
   { /* */
@@ -350,6 +352,8 @@ class ADMaterialProperty : public MaterialPropertyBase<T, true>
 public:
   ADMaterialProperty() = default;
 
+  using typename MaterialPropertyBase<T, true>::value_type;
+
   PropertyValue * init(int size) override { return _init_helper(size, this); }
 
 private:

framework/include/utils/ADOffset.h → framework/include/utils/ADUtils.h

@@ -12,6 +12,12 @@
 #include "MooseError.h"
 #include "MooseTypes.h"
 
+class SystemBase;
+namespace libMesh
+{
+class Elem;
+}
+
 namespace Moose
 {
 /**
@@ -52,33 +58,33 @@ namespace Moose
  * @param max_dofs_per_elem The maximum number of degrees of freedom for any one variable on an
  * element
  * @param element_type The "type" of element that we are on. Current options are
- * \p Moose::ElementType::Element, \p Moose::ElementType::Neighbor, and \p Moose::ElementType::Lower
+ * \p ElementType::Element, \p ElementType::Neighbor, and \p ElementType::Lower
  * @param num_vars_in_system The number of vars in the system. This is used in offset calculation
- * unless \p element_type is \p Moose::ElementType::Element
+ * unless \p element_type is \p ElementType::Element
  * @return The automatic differentiation indexing offset
  *
  */
 inline std::size_t
 adOffset(unsigned int var_num,
          std::size_t max_dofs_per_elem,
-         Moose::ElementType element_type,
+         ElementType element_type = ElementType::Element,
          unsigned int num_vars_in_system = 0)
 {
-  // If our element type is anything other than Moose::ElementType::Element, then the user must
+  // If our element type is anything other than ElementType::Element, then the user must
   // supply num_vars_in_system in order to calculate the offset
-  mooseAssert(element_type == Moose::ElementType::Element || num_vars_in_system,
-              "If our element type is anything other than Moose::ElementType::Element, then you "
+  mooseAssert(element_type == ElementType::Element || num_vars_in_system,
+              "If our element type is anything other than ElementType::Element, then you "
               "must supply num_vars_in_system in order to calculate the offset");
 
   switch (element_type)
   {
-    case Moose::ElementType::Element:
+    case ElementType::Element:
       return var_num * max_dofs_per_elem;
 
-    case Moose::ElementType::Neighbor:
+    case ElementType::Neighbor:
       return num_vars_in_system * max_dofs_per_elem + var_num * max_dofs_per_elem;
 
-    case Moose::ElementType::Lower:
+    case ElementType::Lower:
       return 2 * num_vars_in_system * max_dofs_per_elem + var_num * max_dofs_per_elem;
 
     default:
@@ -91,13 +97,45 @@ adOffset(unsigned int var_num,
 inline std::size_t
 adOffset(unsigned int var_num,
          std::size_t max_dofs_per_elem,
-         Moose::DGJacobianType dg_jacobian_type,
+         DGJacobianType dg_jacobian_type,
          unsigned int num_vars_in_system = 0)
 {
-  if (dg_jacobian_type == Moose::DGJacobianType::ElementElement ||
-      dg_jacobian_type == Moose::DGJacobianType::NeighborElement)
-    return adOffset(var_num, max_dofs_per_elem, Moose::ElementType::Element);
+  if (dg_jacobian_type == DGJacobianType::ElementElement ||
+      dg_jacobian_type == DGJacobianType::NeighborElement)
+    return adOffset(var_num, max_dofs_per_elem, ElementType::Element);
   else
-    return adOffset(var_num, max_dofs_per_elem, Moose::ElementType::Neighbor, num_vars_in_system);
+    return adOffset(var_num, max_dofs_per_elem, ElementType::Neighbor, num_vars_in_system);
+}
+
+/**
+ * Generate a map from global dof index to derivative value
+ */
+std::unordered_map<dof_id_type, Real>
+globalDofIndexToDerivative(const ADReal & ad_real,
+                           const SystemBase & sys,
+                           ElementType elem_type = ElementType::Element,
+                           THREAD_ID tid = 0);
+
+/**
+ * Generate a map from global dof index to derivative value for a (probably quadrature-point-based)
+ * container like a material property or a variable value
+ */
+template <typename T>
+auto
+globalDofIndexToDerivative(const T & ad_real_container,
+                           const SystemBase & sys,
+                           ElementType elem_type = ElementType::Element,
+                           THREAD_ID tid = 0)
+    -> std::vector<std::unordered_map<
+        dof_id_type,
+        typename std::enable_if<std::is_same<ADReal, typename T::value_type>::value, Real>::type>>
+{
+  std::vector<std::unordered_map<dof_id_type, Real>> ret_val(ad_real_container.size());
+
+  for (MooseIndex(ad_real_container) i = 0; i < ad_real_container.size(); ++i)
+    ret_val[i] = globalDofIndexToDerivative(ad_real_container[i], sys, elem_type, tid);
+
+  return ret_val;
 }
+
 }

framework/include/utils/MooseArray.h

@@ -16,6 +16,8 @@ template <typename T>
 class MooseArray
 {
 public:
+  typedef T value_type;
+
   /**
    * Default constructor.  Doesn't initialize anything.
    */

framework/include/variables/MooseVariableData.h

@@ -17,7 +17,7 @@
 #include "libmesh/tensor_value.h"
 #include "libmesh/type_n_tensor.h"
 #include "libmesh/fe_type.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 #include <functional>
 #include <vector>

framework/src/bcs/ADIntegratedBC.C

@@ -15,7 +15,7 @@
 #include "SystemBase.h"
 #include "MooseVariableFE.h"
 #include "MooseVariableScalar.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 #include "libmesh/quadrature.h"
 

framework/src/constraints/ADMortarConstraint.C

@@ -13,7 +13,7 @@
 #include "MooseVariable.h"
 #include "Assembly.h"
 #include "SystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 InputParameters
 ADMortarConstraint::validParams()

framework/src/dgkernels/ADDGKernel.C

@@ -13,7 +13,7 @@
 #include "Problem.h"
 #include "SubProblem.h"
 #include "NonlinearSystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/threads.h"

framework/src/fvbcs/FVFluxBC.C

@@ -11,6 +11,7 @@
 #include "MooseVariableFV.h"
 #include "SystemBase.h"
 #include "Assembly.h"
+#include "ADUtils.h"
 
 InputParameters
 FVFluxBC::validParams()
@@ -85,9 +86,7 @@ FVFluxBC::computeJacobian(Moose::DGJacobianType type, const ADReal & residual)
     SystemBase & sys = _subproblem.systemBaseNonlinear();
     auto dofs_per_elem = sys.getMaxVarNDofsPerElem();
 
-    auto ad_offset = jvar * dofs_per_elem;
-    if (type == Moose::ElementNeighbor || type == Moose::NeighborNeighbor)
-      ad_offset += sys.system().n_vars() * dofs_per_elem;
+    auto ad_offset = Moose::adOffset(jvar, dofs_per_elem, type, sys.system().n_vars());
 
     prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
 

framework/src/fvkernels/FVElementalKernel.C

@@ -12,6 +12,7 @@
 #include "Assembly.h"
 #include "SubProblem.h"
 #include "NonlinearSystemBase.h"
+#include "ADUtils.h"
 
 #include "libmesh/elem.h"
 
@@ -60,7 +61,7 @@ FVElementalKernel::computeJacobian()
 {
   prepareMatrixTag(_assembly, _var.number(), _var.number());
   auto dofs_per_elem = _subproblem.systemBaseNonlinear().getMaxVarNDofsPerElem();
-  size_t ad_offset = _var.number() * dofs_per_elem;
+  auto ad_offset = Moose::adOffset(_var.number(), dofs_per_elem);
   const auto r = computeQpResidual() * _assembly.elemVolume();
   _local_ke(0, 0) += r.derivatives()[ad_offset];
   accumulateTaggedLocalMatrix();
@@ -87,7 +88,8 @@ FVElementalKernel::computeOffDiagJacobian()
     if (ivar != _var.number())
       continue;
 
-    size_t ad_offset = jvar * _subproblem.systemBaseNonlinear().getMaxVarNDofsPerElem();
+    auto ad_offset =
+        Moose::adOffset(jvar, _subproblem.systemBaseNonlinear().getMaxVarNDofsPerElem());
 
     prepareMatrixTag(_assembly, ivar, jvar);
 

framework/src/fvkernels/FVFluxKernel.C

@@ -13,6 +13,7 @@
 #include "SystemBase.h"
 #include "FVDirichletBC.h"
 #include "MooseMesh.h"
+#include "ADUtils.h"
 #include "libmesh/elem.h"
 
 InputParameters
@@ -130,9 +131,7 @@ FVFluxKernel::computeJacobian(Moose::DGJacobianType type, const ADReal & residua
     SystemBase & sys = _subproblem.systemBaseNonlinear();
     auto dofs_per_elem = sys.getMaxVarNDofsPerElem();
 
-    auto ad_offset = jvar * dofs_per_elem;
-    if (type == Moose::ElementNeighbor || type == Moose::NeighborNeighbor)
-      ad_offset += sys.system().n_vars() * dofs_per_elem;
+    auto ad_offset = Moose::adOffset(jvar, dofs_per_elem, type, sys.system().n_vars());
 
     prepareMatrixTagNeighbor(_assembly, ivar, jvar, type);
 

framework/src/interfacekernels/ADInterfaceKernel.C

@@ -13,7 +13,7 @@
 #include "Assembly.h"
 #include "MooseVariableFE.h"
 #include "SystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/quadrature.h"

framework/src/kernels/ADKernel.C

@@ -13,7 +13,7 @@
 #include "Problem.h"
 #include "SubProblem.h"
 #include "NonlinearSystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/threads.h"

framework/src/kernels/ADKernelGrad.C

@@ -11,7 +11,7 @@
 #include "MathUtils.h"
 #include "Assembly.h"
 #include "SystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/threads.h"

framework/src/kernels/ADKernelStabilized.C

@@ -11,6 +11,7 @@
 #include "MathUtils.h"
 #include "Assembly.h"
 #include "SystemBase.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/threads.h"
@@ -68,7 +69,7 @@ ADKernelStabilizedTempl<T>::computeJacobian()
 {
   prepareMatrixTag(_assembly, _var.number(), _var.number());
 
-  size_t ad_offset = _var.number() * _sys.getMaxVarNDofsPerElem();
+  auto ad_offset = Moose::adOffset(_var.number(), _sys.getMaxVarNDofsPerElem());
 
   precalculateResidual();
 
@@ -147,7 +148,7 @@ ADKernelStabilizedTempl<T>::computeADOffDiagJacobian()
     if (ivar != _var.number())
       continue;
 
-    size_t ad_offset = jvar * _sys.getMaxVarNDofsPerElem();
+    auto ad_offset = Moose::adOffset(jvar, _sys.getMaxVarNDofsPerElem());
 
     prepareMatrixTag(_assembly, ivar, jvar);
 

framework/src/kernels/ADKernelValue.C

@@ -10,7 +10,7 @@
 #include "ADKernelValue.h"
 #include "Assembly.h"
 #include "SystemBase.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 // libmesh includes
 #include "libmesh/threads.h"

framework/src/utils/ADUtils.C

@@ -0,0 +1,76 @@
+//* 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 "ADUtils.h"
+
+#include "SystemBase.h"
+#include "SubProblem.h"
+#include "Assembly.h"
+#include "MooseError.h"
+#include "libmesh/system.h"
+#include "libmesh/dof_map.h"
+
+namespace Moose
+{
+
+std::unordered_map<dof_id_type, Real>
+globalDofIndexToDerivative(const ADReal & ad_real,
+                           const SystemBase & sys,
+                           const ElementType elem_type /*=ElementType::Element*/,
+                           const THREAD_ID tid /*=0*/)
+{
+  const Assembly & assembly = sys.subproblem().assembly(tid);
+  const Elem * elem;
+  switch (elem_type)
+  {
+    case ElementType::Element:
+      elem = assembly.elem();
+      break;
+
+    case ElementType::Neighbor:
+      elem = assembly.neighbor();
+      break;
+
+    case ElementType::Lower:
+      elem = assembly.lowerDElem();
+      break;
+
+    default:
+      mooseError("Unrecognized element type");
+  }
+
+  std::unordered_map<dof_id_type, Real> ret_val;
+
+  const System & libmesh_sys = sys.system();
+  const DofMap & dof_map = libmesh_sys.get_dof_map();
+
+  const unsigned int num_vars = libmesh_sys.n_vars();
+
+  const auto max_dofs_per_elem = sys.getMaxVarNDofsPerElem();
+
+  for (unsigned int var_num = 0; var_num < num_vars; ++var_num)
+  {
+    std::vector<dof_id_type> global_indices;
+
+    // Get the global indices corresponding to var_num that exist on elem
+    dof_map.dof_indices(elem, global_indices, var_num);
+
+    // determine the AD offset for the current var
+    const auto ad_offset = adOffset(var_num, max_dofs_per_elem, elem_type, num_vars);
+
+    // Map from global index to derivative
+    for (MooseIndex(global_indices) local_index = 0; local_index < global_indices.size();
+         ++local_index)
+      ret_val[global_indices[local_index]] = ad_real.derivatives()[ad_offset + local_index];
+  }
+
+  return ret_val;
+}
+
+}

framework/src/variables/MooseVariableDataFV.C

@@ -20,7 +20,7 @@
 #include "FVDirichletBC.h"
 #include "SubProblem.h"
 #include "FVKernel.h"
-#include "ADOffset.h"
+#include "ADUtils.h"
 
 #include "libmesh/quadrature.h"
 #include "libmesh/fe_base.h"