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             : #include "PorousFlowMaterial.h"
      11             : 
      12             : #include "MaterialPropertyStorage.h"
      13             : 
      14             : #include "libmesh/quadrature.h"
      15             : 
      16             : #include <limits>
      17             : 
      18             : InputParameters
      19      571434 : PorousFlowMaterial::validParams()
      20             : {
      21      571434 :   InputParameters params = Material::validParams();
      22     1142868 :   params.addRequiredParam<UserObjectName>(
      23             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      24     1142868 :   params.addParam<bool>(
      25     1142868 :       "at_nodes", false, "Evaluate Material properties at nodes instead of quadpoints");
      26     1142868 :   params.addPrivateParam<std::string>("pf_material_type", "pf_material");
      27      571434 :   params.addClassDescription("This generalises MOOSE's Material class to allow for Materials that "
      28             :                              "hold information related to the nodes in the finite element");
      29             : 
      30             :   // Needed due to the custom tomfoolery going on with nodal material sizing in
      31             :   // initStatefulProperties()
      32      571434 :   params.set<bool>("_force_stateful_init") = true;
      33             : 
      34      571434 :   return params;
      35           0 : }
      36             : 
      37      447095 : PorousFlowMaterial::PorousFlowMaterial(const InputParameters & parameters)
      38             :   : Material(parameters),
      39      447095 :     _nodal_material(getParam<bool>("at_nodes")),
      40      447095 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      41             :     _pressure_variable_name("pressure_variable"),
      42             :     _saturation_variable_name("saturation_variable"),
      43             :     _temperature_variable_name("temperature_variable"),
      44      447095 :     _mass_fraction_variable_name("mass_fraction_variable")
      45             : {
      46      447095 : }
      47             : 
      48             : void
      49      433821 : PorousFlowMaterial::initialSetup()
      50             : {
      51      433821 :   if (!_nodal_material)
      52             :     return;
      53             : 
      54      197829 :   _material_data.onlyResizeIfSmaller(true);
      55             :   auto & storage = _material_data.getMaterialPropertyStorage();
      56      197829 :   if (!storage.hasStatefulProperties())
      57             :     return;
      58             : 
      59             :   auto & stateful_prop_id_to_prop_id = storage.statefulProps();
      60      533607 :   for (const auto i : index_range(stateful_prop_id_to_prop_id))
      61             :   {
      62      438591 :     const auto prop_id = stateful_prop_id_to_prop_id[i];
      63             :     if (_supplied_prop_ids.count(prop_id))
      64       29692 :       _supplied_old_prop_ids.push_back(i);
      65             :   }
      66             : }
      67             : 
      68             : void
      69     7093661 : PorousFlowMaterial::initStatefulProperties(unsigned int n_points)
      70             : {
      71     7093661 :   if (_nodal_material)
      72             :   {
      73             :     // size the properties to max(number_of_nodes, number_of_quadpoints)
      74     3131368 :     sizeNodalProperties();
      75             : 
      76             :     // compute the values for number_of_nodes
      77     3131368 :     Material::initStatefulProperties(_current_elem->n_nodes());
      78             :   }
      79             :   else
      80     3962293 :     Material::initStatefulProperties(n_points);
      81     7093661 : }
      82             : 
      83             : void
      84   106171382 : PorousFlowMaterial::computeNodalProperties()
      85             : {
      86   106171382 :   const unsigned int numnodes = _current_elem->n_nodes();
      87             : 
      88             :   // compute the values for all nodes
      89   534270630 :   for (_qp = 0; _qp < numnodes; ++_qp)
      90   428099330 :     computeQpProperties();
      91             : 
      92             :   // If number_of_nodes < number_of_quadpoints, the remaining values in the
      93             :   // material data array are zero (for scalars) and empty (for vectors).
      94             :   // Unfortunately, this can cause issues with adaptivity, where the empty
      95             :   // value can be transferred to a node in a child element. This can lead
      96             :   // to a segfault when accessing stateful properties, see #14428.
      97             :   // To prevent this, we copy the last node value to the empty array positions.
      98   106171300 :   if (numnodes < _qrule->n_points())
      99             :   {
     100    11766624 :     MaterialProperties & props = _material_data.props();
     101             : 
     102             :     // Copy from qp = _current_elem->n_nodes() - 1 to qp = _qrule->n_points() -1
     103    46478988 :     for (const auto & prop_id : _supplied_prop_ids)
     104   136574912 :       for (unsigned int qp = numnodes; qp < _qrule->n_points(); ++qp)
     105   101862548 :         props[prop_id].qpCopy(qp, props[prop_id], numnodes - 1);
     106             :   }
     107   106171300 : }
     108             : 
     109             : void
     110   218912318 : PorousFlowMaterial::computeProperties()
     111             : {
     112   218912318 :   if (_nodal_material)
     113             :   {
     114             :     // size the properties to max(number_of_nodes, number_of_quadpoints)
     115   106171382 :     sizeNodalProperties();
     116             : 
     117   106171382 :     computeNodalProperties();
     118             :   }
     119             :   else
     120   112740936 :     Material::computeProperties();
     121   218912182 : }
     122             : 
     123             : void
     124   109302750 : PorousFlowMaterial::sizeNodalProperties()
     125             : {
     126             :   /*
     127             :    * For nodal materials, the Properties should be sized as the maximum of
     128             :    * the number of nodes and the number of quadpoints.
     129             :    * We only actually need "number of nodes" pieces of information, which are
     130             :    * computed by computeProperties(), so the n_points - _current_elem->n_nodes()
     131             :    * elements at the end of the std::vector will always be zero, but they
     132             :    * are needed because MOOSE does copy operations (etc) that assumes that
     133             :    * the std::vector is sized to number of quadpoints.
     134             :    *
     135             :    * On boundary materials, the number of nodes may be larger than the number of
     136             :    * qps on the face of the element, in which case the remaining entries in the
     137             :    * material properties storage will be zero.
     138             :    *
     139             :    * \author lindsayad: MooseArray currently has the unfortunate side effect that if your new size
     140             :    * is greater than the current size, then we clear the whole data structure. Consequently this
     141             :    * call has the potential to clear material property evaluations done earlier in the material
     142             :    * dependency chain. So instead we selectively resize just our own properties and not everyone's
     143             :    */
     144             :   // _material_data.resize(std::max(_current_elem->n_nodes(), _qrule->n_points()));
     145             : 
     146   109302750 :   const auto new_size = std::max(_current_elem->n_nodes(), _qrule->n_points());
     147   109302750 :   auto & storage = _material_data.getMaterialPropertyStorage();
     148             : 
     149             :   auto & props = _material_data.props();
     150   454690225 :   for (const auto prop_id : _supplied_prop_ids)
     151   345387475 :     props[prop_id].resize(new_size);
     152             : 
     153   211444706 :   for (const auto state : storage.statefulIndexRange())
     154   146193752 :     for (const auto prop_id : _supplied_old_prop_ids)
     155    44051796 :       if (_material_data.props(state).hasValue(prop_id))
     156     2157085 :         _material_data.props(state)[prop_id].resize(new_size);
     157   109302750 : }
     158             : 
     159             : unsigned
     160      341062 : PorousFlowMaterial::nearestQP(unsigned nodenum) const
     161             : {
     162             :   unsigned nearest_qp = 0;
     163             :   Real smallest_dist = std::numeric_limits<Real>::max();
     164     1875994 :   for (const auto qp : make_range(_qrule->n_points()))
     165             :   {
     166     1534932 :     const Real this_dist = (_current_elem->point(nodenum) - _q_point[qp]).norm();
     167     1534932 :     if (this_dist < smallest_dist)
     168             :     {
     169             :       nearest_qp = qp;
     170             :       smallest_dist = this_dist;
     171             :     }
     172             :   }
     173      341062 :   return nearest_qp;
     174             : }