Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "MooseError.h"
      13             : #include "MooseTypes.h"
      14             : #include "MooseConfig.h"
      15             : 
      16             : class SubProblem;
      17             : class SystemBase;
      18             : namespace libMesh
      19             : {
      20             : class Elem;
      21             : }
      22             : 
      23             : namespace Moose
      24             : {
      25             : /**
      26             :  * Whether we are using global AD indexing
      27             :  */
      28             : inline bool
      29        6824 : globalADIndexing()
      30             : {
      31        6824 :   return true;
      32             : }
      33             : 
      34             : /**
      35             :  * Helper function for computing automatic differentiation offset. Let's explain how our derivative
      36             :  * index numbering scheme works:
      37             :  *
      38             :  * Let's just think about continuous finite elements for a second. We use a variable major
      39             :  * numbering scheme, such that each variables indices are in a contiguous block. Let's imagine we
      40             :  * have two variables, u and v, and we're on a \p QUAD4. The AD indices will be ordered like this:
      41             :  *
      42             :  * u0, u1, u2, u3, v0, v1, v2, v3
      43             :  *
      44             :  * \p max_dofs_per_elem should be for a \p QUAD4: 4. For a \p QUAD9, 9. \p HEX27, 27. Etc. For
      45             :  * CFEM the offset will be simply be the \p max_dofs_per_elem number times the \p var_num. So for u
      46             :  * on a \p QUAD4: 4 * 0 = 0. For v: 4 * 1. So u indices start at index 0, v indices start at
      47             :  * index 4.
      48             :  *
      49             :  * With DFEM or interface kernels it's a little more complicated. We essentially already have an
      50             :  * indices block that is \p num_vars_in_system * \p max_dofs_per_elem long, so in our two var,
      51             :  * \p QUAD4 example: 4 * 2 = 8. So what we do is that if we are on a neighbor element, we do an
      52             :  * additional offset by \p num_vars_in_system * \p max_dofs_per_elem. So now our derivative indices
      53             :  * are ordered like this:
      54             :  *
      55             :  * u0, u1, u2, u3, v0, v1, v2, v3, u0_neighbor, u1_neighbor, u2_neighbor, u3_neighbor, v0_neighbor,
      56             :  * v1_neighbor, v2_neighbor, v3_neighbor
      57             :  *
      58             :  * Finally if a lower-dimensional element is involved, then we another offset of \p
      59             :  * num_vars_in_system * \p max_dofs_per_elem:
      60             :  *
      61             :  * u0, u1, u2, u3, v0, v1, v2, v3, u0_neighbor, u1_neighbor, u2_neighbor, u3_neighbor, v0_neighbor,
      62             :  * v1_neighbor, v2_neighbor, v3_neighbor, u0_lower, u1_lower, u2_lower, u3_lower, v0_lower,
      63             :  * v1_lower, v2_lower, v3_lower
      64             :  *
      65             :  * Note that a lower dimensional block will have less indices than a higher dimensional one, but we
      66             :  * do not optimize for that consideration at this time
      67             :  *
      68             :  * @param var_num The variable number we are calculating the offset for
      69             :  * @param max_dofs_per_elem The maximum number of degrees of freedom for any one variable on an
      70             :  * element
      71             :  * @param element_type The "type" of element that we are on. Current options are
      72             :  * \p ElementType::Element, \p ElementType::Neighbor, and \p ElementType::Lower
      73             :  * @param num_vars_in_system The number of vars in the system. This is used in offset calculation
      74             :  * unless \p element_type is \p ElementType::Element
      75             :  * @return The automatic differentiation indexing offset
      76             :  *
      77             :  */
      78             : inline std::size_t
      79           0 : adOffset(unsigned int var_num,
      80             :          std::size_t max_dofs_per_elem,
      81             :          ElementType element_type = ElementType::Element,
      82             :          unsigned int num_vars_in_system = 0)
      83             : {
      84             :   // If our element type is anything other than ElementType::Element, then the user must
      85             :   // supply num_vars_in_system in order to calculate the offset
      86             :   mooseAssert(element_type == ElementType::Element || num_vars_in_system,
      87             :               "If our element type is anything other than ElementType::Element, then you "
      88             :               "must supply num_vars_in_system in order to calculate the offset");
      89             : 
      90           0 :   switch (element_type)
      91             :   {
      92           0 :     case ElementType::Element:
      93           0 :       return var_num * max_dofs_per_elem;
      94             : 
      95           0 :     case ElementType::Neighbor:
      96           0 :       return num_vars_in_system * max_dofs_per_elem + var_num * max_dofs_per_elem;
      97             : 
      98           0 :     case ElementType::Lower:
      99           0 :       return 2 * num_vars_in_system * max_dofs_per_elem + var_num * max_dofs_per_elem;
     100             : 
     101           0 :     default:
     102           0 :       mooseError(
     103             :           "Unsupported element type ",
     104           0 :           static_cast<typename std::underlying_type<decltype(element_type)>::type>(element_type));
     105             :   }
     106             : }
     107             : 
     108             : inline std::size_t
     109             : adOffset(unsigned int var_num,
     110             :          std::size_t max_dofs_per_elem,
     111             :          DGJacobianType dg_jacobian_type,
     112             :          unsigned int num_vars_in_system = 0)
     113             : {
     114             :   if (dg_jacobian_type == DGJacobianType::ElementElement ||
     115             :       dg_jacobian_type == DGJacobianType::NeighborElement)
     116             :     return adOffset(var_num, max_dofs_per_elem, ElementType::Element);
     117             :   else
     118             :     return adOffset(var_num, max_dofs_per_elem, ElementType::Neighbor, num_vars_in_system);
     119             : }
     120             : 
     121             : /**
     122             :  * Generate a map from global dof index to derivative value
     123             :  */
     124             : std::unordered_map<dof_id_type, Real>
     125             : globalDofIndexToDerivative(const ADReal & ad_real,
     126             :                            const SystemBase & sys,
     127             :                            ElementType elem_type = ElementType::Element,
     128             :                            THREAD_ID tid = 0);
     129             : 
     130             : /**
     131             :  * Generate a map from global dof index to derivative value for a (probably quadrature-point-based)
     132             :  * container like a material property or a variable value
     133             :  */
     134             : template <typename T>
     135             : auto
     136             : globalDofIndexToDerivative(const T & ad_real_container,
     137             :                            const SystemBase & sys,
     138             :                            ElementType elem_type = ElementType::Element,
     139             :                            THREAD_ID tid = 0)
     140             :     -> std::vector<std::unordered_map<
     141             :         dof_id_type,
     142             :         typename std::enable_if<std::is_same<ADReal, typename T::value_type>::value, Real>::type>>
     143             : {
     144             :   std::vector<std::unordered_map<dof_id_type, Real>> ret_val(ad_real_container.size());
     145             : 
     146             :   for (MooseIndex(ad_real_container) i = 0; i < ad_real_container.size(); ++i)
     147             :     ret_val[i] = globalDofIndexToDerivative(ad_real_container[i], sys, elem_type, tid);
     148             : 
     149             :   return ret_val;
     150             : }
     151             : 
     152             : /**
     153             :  * @returns whether we should be doing derivatives
     154             :  */
     155             : bool doDerivatives(const SubProblem & subproblem, const SystemBase & sys);
     156             : }