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 "MooseArray.h"
      13             : #include "MooseTypes.h"
      14             : #include "MeshChangedInterface.h"
      15             : #include "MooseVariableDataBase.h"
      16             : #include "TheWarehouse.h"
      17             : 
      18             : #include "libmesh/tensor_tools.h"
      19             : #include "libmesh/vector_value.h"
      20             : #include "libmesh/tensor_value.h"
      21             : #include "libmesh/type_n_tensor.h"
      22             : #include "libmesh/fe_type.h"
      23             : #include "libmesh/dof_map.h"
      24             : #include "libmesh/enum_fe_family.h"
      25             : #include "SubProblem.h"
      26             : 
      27             : #include <functional>
      28             : #include <vector>
      29             : 
      30             : class FaceInfo;
      31             : class SystemBase;
      32             : class TimeIntegrator;
      33             : class Assembly;
      34             : 
      35             : template <typename>
      36             : class MooseVariableFV;
      37             : 
      38             : namespace libMesh
      39             : {
      40             : class QBase;
      41             : }
      42             : 
      43             : namespace Moose
      44             : {
      45             : template <typename T>
      46             : void
      47   895177888 : initDofIndices(T & data, const Elem & elem)
      48             : {
      49   895177888 :   if (data._prev_elem != &elem)
      50             :   {
      51   226541830 :     data._dof_map.dof_indices(&elem, data._dof_indices, data._var_num);
      52   226541830 :     data._prev_elem = &elem;
      53             :   }
      54   895177888 : }
      55             : }
      56             : 
      57             : namespace
      58             : {
      59             : template <typename T, typename T2>
      60             : void
      61   125897451 : assignForAllQps(const T & value, T2 & array, const unsigned int nqp)
      62             : {
      63   253925341 :   for (const auto qp : make_range(nqp))
      64   128027890 :     array[qp] = value;
      65   125897451 : }
      66             : }
      67             : 
      68             : template <typename OutputType>
      69             : class MooseVariableDataFV : public MooseVariableDataBase<OutputType>, public MeshChangedInterface
      70             : {
      71             : public:
      72             :   // type for gradient, second and divergence of template class OutputType
      73             :   typedef typename libMesh::TensorTools::IncrementRank<OutputType>::type OutputGradient;
      74             :   typedef typename libMesh::TensorTools::IncrementRank<OutputGradient>::type OutputSecond;
      75             :   typedef typename libMesh::TensorTools::DecrementRank<OutputType>::type OutputDivergence;
      76             : 
      77             :   // shortcut for types storing values on quadrature points
      78             :   typedef MooseArray<OutputType> FieldVariableValue;
      79             :   typedef MooseArray<OutputGradient> FieldVariableGradient;
      80             :   typedef MooseArray<OutputSecond> FieldVariableSecond;
      81             :   typedef MooseArray<OutputType> FieldVariableCurl;
      82             :   typedef MooseArray<OutputDivergence> FieldVariableDivergence;
      83             : 
      84             :   // shape function type for the template class OutputType
      85             :   typedef typename Moose::ShapeType<OutputType>::type OutputShape;
      86             : 
      87             :   // type for gradient, second and divergence of shape functions of template class OutputType
      88             :   typedef typename libMesh::TensorTools::IncrementRank<OutputShape>::type OutputShapeGradient;
      89             :   typedef typename libMesh::TensorTools::IncrementRank<OutputShapeGradient>::type OutputShapeSecond;
      90             :   typedef typename libMesh::TensorTools::DecrementRank<OutputShape>::type OutputShapeDivergence;
      91             : 
      92             :   // DoF value type for the template class OutputType
      93             :   typedef typename Moose::DOFType<OutputType>::type OutputData;
      94             :   typedef MooseArray<OutputData> DofValue;
      95             : 
      96             :   MooseVariableDataFV(const MooseVariableFV<OutputType> & var,
      97             :                       SystemBase & sys,
      98             :                       THREAD_ID tid,
      99             :                       Moose::ElementType element_type,
     100             :                       const Elem * const & elem);
     101             : 
     102           0 :   bool isNodal() const override { return false; }
     103           0 :   bool hasDoFsOnNodes() const override { return false; }
     104         344 :   libMesh::FEContinuity getContinuity() const override { return libMesh::DISCONTINUOUS; }
     105             : 
     106             :   /**
     107             :    * Returns whether this data structure needs automatic differentiation calculations
     108             :    */
     109           0 :   bool needsAD() const { return _need_ad; }
     110             : 
     111             :   /**
     112             :    * Set the geometry type before calculating variables values
     113             :    * @param gm_type The type type of geometry; either Volume or Face
     114             :    */
     115             :   void setGeometry(Moose::GeometryType gm_type);
     116             : 
     117             :   //////////////// Heavy lifting computational routines //////////////////////////////
     118             : 
     119             :   /**
     120             :    * compute the variable values
     121             :    */
     122             :   void computeValuesFace(const FaceInfo & fi);
     123             : 
     124             :   /**
     125             :    * compute the variable values
     126             :    */
     127             :   void computeValues();
     128             : 
     129             :   /**
     130             :    * compute AD things
     131             :    */
     132             :   void computeAD(const unsigned int num_dofs, const unsigned int nqp);
     133             : 
     134             :   ///////////////////////////// Shape functions /////////////////////////////////////
     135             : 
     136             :   /**
     137             :    * The current element
     138             :    */
     139         664 :   const Elem * const & currentElem() const { return _elem; }
     140             : 
     141             :   /**
     142             :    * prepare the initial condition
     143             :    */
     144             :   void prepareIC();
     145             : 
     146             :   //////////////////////////////////// Solution getters /////////////////////////////////////
     147             : 
     148             :   /**
     149             :    * Local solution value
     150             :    */
     151             :   const FieldVariableValue & sln(Moose::SolutionState state) const;
     152             : 
     153             :   /**
     154             :    * Local time derivative of solution gradient getter
     155             :    */
     156             :   const FieldVariableGradient & gradSlnDot() const;
     157             : 
     158             :   /**
     159             :    * Local second time derivative of solution gradient getter
     160             :    */
     161             :   const FieldVariableGradient & gradSlnDotDot() const;
     162             : 
     163             :   /**
     164             :    * Local solution second spatial derivative getter
     165             :    * @param state The state of the simulation: current, old, older, previous nl
     166             :    */
     167             :   const FieldVariableSecond & secondSln(Moose::SolutionState state) const;
     168             : 
     169             :   /**
     170             :    * Local solution curl getter
     171             :    * @param state The state of the simulation: current, old, older
     172             :    */
     173             :   const FieldVariableCurl & curlSln(Moose::SolutionState state) const;
     174             : 
     175             :   /**
     176             :    * @returns The automatic differentiation solution indexable at quadrature points
     177             :    */
     178             :   const ADTemplateVariableValue<OutputType> & adSln() const;
     179             : 
     180             :   /**
     181             :    * @returns The automatic differentiation gradient indexable at quadrature points
     182             :    */
     183             :   const ADTemplateVariableGradient<OutputType> & adGradSln() const;
     184             : 
     185           0 :   const ADTemplateVariableGradient<OutputType> & adGradSlnDot() const
     186             :   {
     187           0 :     mooseError("Gradient of time derivative not yet implemented for FV");
     188             :   }
     189             : 
     190             :   /**
     191             :    * @returns The automatic differentiation matrix of second derivatives indexable at quadrature
     192             :    * points
     193             :    */
     194             :   const ADTemplateVariableSecond<OutputType> & adSecondSln() const;
     195             : 
     196             :   const ADTemplateVariableValue<OutputType> & adUDot() const;
     197             : 
     198             :   const ADTemplateVariableValue<OutputType> & adUDotDot() const;
     199             : 
     200             :   const FieldVariableValue & uDot() const;
     201             : 
     202             :   const FieldVariableValue & uDotDot() const;
     203             : 
     204             :   const FieldVariableValue & uDotOld() const;
     205             : 
     206             :   const FieldVariableValue & uDotDotOld() const;
     207             : 
     208           0 :   const VariableValue & duDotDu() const
     209             :   {
     210           0 :     _need_du_dot_du = true;
     211           0 :     return _du_dot_du;
     212             :   }
     213             : 
     214           0 :   const VariableValue & duDotDotDu() const
     215             :   {
     216           0 :     _need_du_dotdot_du = true;
     217           0 :     return _du_dotdot_du;
     218             :   }
     219             : 
     220             :   /**
     221             :    * Set local DOF values and evaluate the values on quadrature points
     222             :    */
     223             :   void setDofValues(const DenseVector<OutputData> & values);
     224             : 
     225             :   ///@{
     226             :   /**
     227             :    * dof value setters
     228             :    */
     229             :   void setDofValue(const OutputData & value, unsigned int index);
     230             :   ///@}
     231             : 
     232             :   OutputData
     233             :   getElementalValue(const Elem * elem, Moose::SolutionState state, unsigned int idx = 0) const;
     234             : 
     235             :   ///////////////////////////// dof indices ///////////////////////////////////////////////
     236             : 
     237             :   void getDofIndices(const Elem * elem, std::vector<dof_id_type> & dof_indices) const;
     238             :   const std::vector<dof_id_type> & dofIndices() const;
     239             :   unsigned int numberOfDofs() const;
     240    10686570 :   void clearDofIndices()
     241             :   {
     242    10686570 :     _dof_indices.clear();
     243    10686570 :     _prev_elem = nullptr;
     244    10686570 :   }
     245             : 
     246             :   /////////////////////////// DoF value getters /////////////////////////////////////
     247             : 
     248             :   const DofValue & dofValuesDot() const;
     249             :   const DofValue & dofValuesDotOld() const;
     250             :   const DofValue & dofValuesDotDot() const;
     251             :   const DofValue & dofValuesDotDotOld() const;
     252             :   const MooseArray<libMesh::Number> & dofValuesDuDotDu() const;
     253             :   const MooseArray<libMesh::Number> & dofValuesDuDotDotDu() const;
     254             : 
     255             :   /**
     256             :    * Return the AD dof values
     257             :    */
     258             :   const MooseArray<ADReal> & adDofValues() const;
     259             : 
     260             :   /**
     261             :    * Return the AD dof time derivatives
     262             :    */
     263             :   const MooseArray<ADReal> & adDofValuesDot() const;
     264             : 
     265             :   /////////////////////////////// Increment stuff ///////////////////////////////////////
     266             : 
     267             :   /**
     268             :    * Increment getter
     269             :    * @return The increment
     270             :    */
     271           0 :   const FieldVariableValue & increment() const { return _increment; }
     272             : 
     273             :   /**
     274             :    * Compute and store incremental change in solution at QPs based on increment_vec
     275             :    */
     276             :   void computeIncrementAtQps(const libMesh::NumericVector<libMesh::Number> & increment_vec);
     277             : 
     278             :   /// checks if a Dirichlet BC exists on this face
     279           0 :   bool hasDirichletBC() const { return _has_dirichlet_bc; }
     280             : 
     281             :   void meshChanged() override;
     282             : 
     283             : protected:
     284           0 :   virtual const MooseVariableFV<OutputType> & var() const override { return _var; }
     285             : 
     286             : private:
     287             :   void initializeSolnVars();
     288             : 
     289             :   /**
     290             :    * Helper methods for assigning nodal values from their corresponding solution values (dof
     291             :    * values as they're referred to here in this class). These methods are only truly meaningful
     292             :    * for nodal basis families
     293             :    */
     294             :   void fetchADDofValues();
     295             : 
     296             :   /**
     297             :    * Helper method that tells us whether it's safe to compute _ad_u_dot
     298             :    */
     299             :   bool safeToComputeADUDot() const;
     300             : 
     301             :   /// A const reference to the owning MooseVariableFV object
     302             :   const MooseVariableFV<OutputType> & _var;
     303             : 
     304             :   const libMesh::FEType & _fe_type;
     305             : 
     306             :   const unsigned int _var_num;
     307             : 
     308             :   const Assembly & _assembly;
     309             : 
     310             :   /// The element type this object is storing data for. This is either Element, Neighbor, or Lower
     311             :   Moose::ElementType _element_type;
     312             : 
     313             :   /// Continuity type of the variable
     314             :   libMesh::FEContinuity _continuity;
     315             : 
     316             :   /// Increment in the variable used in dampers
     317             :   FieldVariableValue _increment;
     318             : 
     319             :   /// A zero AD variable
     320             :   const ADReal _ad_zero;
     321             : 
     322             :   /// SolutionState second_u flags
     323             :   mutable bool _need_second;
     324             :   mutable bool _need_second_old;
     325             :   mutable bool _need_second_older;
     326             :   mutable bool _need_second_previous_nl;
     327             : 
     328             :   /// curl flags
     329             :   mutable bool _need_curl;
     330             :   mutable bool _need_curl_old;
     331             :   mutable bool _need_curl_older;
     332             : 
     333             :   /// AD flags
     334             :   mutable bool _need_ad;
     335             :   mutable bool _need_ad_u;
     336             :   mutable bool _need_ad_u_dot;
     337             :   mutable bool _need_ad_u_dotdot;
     338             :   mutable bool _need_ad_grad_u;
     339             :   mutable bool _need_ad_grad_u_dot;
     340             :   mutable bool _need_ad_second_u;
     341             : 
     342             :   /// grad_u dots
     343             :   FieldVariableGradient _grad_u_dot;
     344             :   FieldVariableGradient _grad_u_dotdot;
     345             : 
     346             :   /// second_u
     347             :   FieldVariableSecond _second_u;
     348             :   FieldVariableSecond _second_u_old;
     349             :   FieldVariableSecond _second_u_older;
     350             :   FieldVariableSecond _second_u_previous_nl;
     351             : 
     352             :   /// curl_u
     353             :   FieldVariableCurl _curl_u;
     354             :   FieldVariableCurl _curl_u_old;
     355             :   FieldVariableCurl _curl_u_older;
     356             : 
     357             :   /// AD u
     358             :   ADTemplateVariableValue<OutputShape> _ad_u;
     359             :   ADTemplateVariableGradient<OutputShape> _ad_grad_u;
     360             :   ADTemplateVariableSecond<OutputShape> _ad_second_u;
     361             :   MooseArray<ADReal> _ad_dof_values;
     362             :   MooseArray<ADReal> _ad_dofs_dot;
     363             :   MooseArray<ADReal> _ad_dofs_dotdot;
     364             :   ADTemplateVariableValue<OutputShape> _ad_u_dot;
     365             :   ADTemplateVariableValue<OutputShape> _ad_u_dotdot;
     366             :   ADTemplateVariableGradient<OutputShape> _ad_grad_u_dot;
     367             : 
     368             :   // time derivatives
     369             : 
     370             :   /// u_dot (time derivative)
     371             :   FieldVariableValue _u_dot;
     372             : 
     373             :   /// u_dotdot (second time derivative)
     374             :   FieldVariableValue _u_dotdot;
     375             : 
     376             :   /// u_dot_old (time derivative)
     377             :   FieldVariableValue _u_dot_old;
     378             : 
     379             :   /// u_dotdot_old (second time derivative)
     380             :   FieldVariableValue _u_dotdot_old;
     381             : 
     382             :   /// derivative of u_dot wrt u
     383             :   VariableValue _du_dot_du;
     384             : 
     385             :   /// derivative of u_dotdot wrt u
     386             :   VariableValue _du_dotdot_du;
     387             : 
     388             :   /// Pointer to time integrator
     389             :   const TimeIntegrator * const _time_integrator;
     390             : 
     391             :   /// The current elem. This has to be a reference because the current elem will be constantly
     392             :   /// changing. If we initialized this to point to one elem, then in the next calculation we would
     393             :   /// be pointing to the wrong place!
     394             :   const Elem * const & _elem;
     395             :   /// used to keep track of when dof indices are out of date
     396             :   mutable const Elem * _prev_elem = nullptr;
     397             : 
     398             :   const std::vector<dof_id_type> & initDofIndices();
     399             : 
     400             :   /// if this variable has a dirichlet bc defined on a particular face
     401             :   bool _has_dirichlet_bc;
     402             : 
     403             :   /// Whether this variable is being calculated on a displaced system
     404             :   const bool _displaced;
     405             : 
     406             :   /// The quadrature rule
     407             :   const libMesh::QBase * _qrule;
     408             : 
     409             :   /// A dummy ADReal variable
     410       10854 :   ADReal _ad_real_dummy = 0;
     411             : 
     412             :   /// Cached warehouse query for FVElementalKernels
     413             :   TheWarehouse::QueryCache<> _fv_elemental_kernel_query_cache;
     414             :   /// Cached warehouse query for FVFluxKernels
     415             :   TheWarehouse::QueryCache<> _fv_flux_kernel_query_cache;
     416             : 
     417             :   using MooseVariableDataBase<OutputType>::_sys;
     418             :   using MooseVariableDataBase<OutputType>::_subproblem;
     419             :   using MooseVariableDataBase<OutputType>::_need_vector_tag_dof_u;
     420             :   using MooseVariableDataBase<OutputType>::_need_matrix_tag_dof_u;
     421             :   using MooseVariableDataBase<OutputType>::_vector_tags_dof_u;
     422             :   using MooseVariableDataBase<OutputType>::_matrix_tags_dof_u;
     423             :   using MooseVariableDataBase<OutputType>::_vector_tag_u;
     424             :   using MooseVariableDataBase<OutputType>::_need_vector_tag_u;
     425             :   using MooseVariableDataBase<OutputType>::_vector_tag_grad;
     426             :   using MooseVariableDataBase<OutputType>::_need_vector_tag_grad;
     427             :   using MooseVariableDataBase<OutputType>::_matrix_tag_u;
     428             :   using MooseVariableDataBase<OutputType>::_need_matrix_tag_u;
     429             :   using MooseVariableDataBase<OutputType>::_dof_indices;
     430             :   using MooseVariableDataBase<OutputType>::_has_dof_values;
     431             :   using MooseVariableDataBase<OutputType>::fetchDofValues;
     432             :   using MooseVariableDataBase<OutputType>::assignNodalValue;
     433             :   using MooseVariableDataBase<OutputType>::zeroSizeDofValues;
     434             :   using MooseVariableDataBase<OutputType>::_solution_tag;
     435             :   using MooseVariableDataBase<OutputType>::_old_solution_tag;
     436             :   using MooseVariableDataBase<OutputType>::_older_solution_tag;
     437             :   using MooseVariableDataBase<OutputType>::_previous_nl_solution_tag;
     438             :   using MooseVariableDataBase<OutputType>::_dof_map;
     439             :   using MooseVariableDataBase<OutputType>::_need_u_dot;
     440             :   using MooseVariableDataBase<OutputType>::_need_u_dotdot;
     441             :   using MooseVariableDataBase<OutputType>::_need_u_dot_old;
     442             :   using MooseVariableDataBase<OutputType>::_need_u_dotdot_old;
     443             :   using MooseVariableDataBase<OutputType>::_need_du_dot_du;
     444             :   using MooseVariableDataBase<OutputType>::_need_du_dotdot_du;
     445             :   using MooseVariableDataBase<OutputType>::_need_grad_dot;
     446             :   using MooseVariableDataBase<OutputType>::_need_grad_dotdot;
     447             :   using MooseVariableDataBase<OutputType>::_need_dof_values_dot;
     448             :   using MooseVariableDataBase<OutputType>::_need_dof_values_dotdot;
     449             :   using MooseVariableDataBase<OutputType>::_need_dof_values_dot_old;
     450             :   using MooseVariableDataBase<OutputType>::_need_dof_values_dotdot_old;
     451             :   using MooseVariableDataBase<OutputType>::_need_dof_du_dot_du;
     452             :   using MooseVariableDataBase<OutputType>::_need_dof_du_dotdot_du;
     453             :   using MooseVariableDataBase<OutputType>::_dof_values_dot;
     454             :   using MooseVariableDataBase<OutputType>::_dof_values_dotdot;
     455             :   using MooseVariableDataBase<OutputType>::_dof_values_dot_old;
     456             :   using MooseVariableDataBase<OutputType>::_dof_values_dotdot_old;
     457             :   using MooseVariableDataBase<OutputType>::_dof_du_dot_du;
     458             :   using MooseVariableDataBase<OutputType>::_dof_du_dotdot_du;
     459             :   using MooseVariableDataBase<OutputType>::_tid;
     460             :   using MooseVariableDataBase<OutputType>::_nodal_value_dot;
     461             :   using MooseVariableDataBase<OutputType>::_nodal_value_dotdot;
     462             :   using MooseVariableDataBase<OutputType>::_nodal_value_dot_old;
     463             :   using MooseVariableDataBase<OutputType>::_nodal_value_dotdot_old;
     464             :   using MooseVariableDataBase<OutputType>::_required_vector_tags;
     465             : 
     466             :   friend void Moose::initDofIndices<>(MooseVariableDataFV<OutputType> &, const Elem &);
     467             : };
     468             : 
     469             : /////////////////////// General template definitions //////////////////////////////////////
     470             : 
     471             : template <typename OutputType>
     472             : inline bool
     473    29085924 : MooseVariableDataFV<OutputType>::safeToComputeADUDot() const
     474             : {
     475             :   // If we don't have a time integrator then we have no way to calculate _ad_u_dot because we rely
     476             :   // on calls to TimeIntegrator::computeADTimeDerivatives. Another potential situation where
     477             :   // _ad_u_dot computation is potentially troublesome is if we are an auxiliary variable which uses
     478             :   // the auxiliary system copy of the time integrator. Some derived time integrator classes do setup
     479             :   // in their solve() method, and that solve() method only happens for the nonlinear system copy of
     480             :   // the time integrator.
     481    29085924 :   return _time_integrator && (_var.kind() == Moose::VAR_SOLVER);
     482             : }
     483             : 
     484             : template <typename OutputType>
     485             : const ADTemplateVariableValue<OutputType> &
     486          13 : MooseVariableDataFV<OutputType>::adUDotDot() const
     487             : {
     488             :   // Generally speaking, we need u dot information when computing u dot dot
     489          13 :   adUDot();
     490             : 
     491          13 :   _need_ad = _need_ad_u_dotdot = true;
     492             : 
     493          13 :   if (!safeToComputeADUDot())
     494             :     // We will just copy the value of _u_dotdot into _ad_u_dotdot
     495           0 :     _need_u_dotdot = true;
     496             : 
     497          13 :   return _ad_u_dotdot;
     498             : }
     499             : 
     500             : template <typename OutputType>
     501             : const std::vector<dof_id_type> &
     502   640215530 : MooseVariableDataFV<OutputType>::dofIndices() const
     503             : {
     504   640215530 :   return const_cast<MooseVariableDataFV<OutputType> *>(this)->initDofIndices();
     505             : }
     506             : 
     507             : template <typename OutputType>
     508             : unsigned int
     509      177059 : MooseVariableDataFV<OutputType>::numberOfDofs() const
     510             : {
     511      177059 :   return dofIndices().size();
     512             : }