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