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 "MeshDivisionFunctorReductionVectorPostprocessor.h"
      11             : #include "MeshDivision.h"
      12             : #include "MooseMesh.h"
      13             : #include "MooseMeshUtils.h"
      14             : #include "FEProblemBase.h"
      15             : 
      16             : #include "libmesh/mesh_base.h"
      17             : 
      18             : registerMooseObject("MooseApp", MeshDivisionFunctorReductionVectorPostprocessor);
      19             : 
      20             : InputParameters
      21        3245 : MeshDivisionFunctorReductionVectorPostprocessor::validParams()
      22             : {
      23        3245 :   InputParameters params = SpatialUserObjectFunctor<ElementVectorPostprocessor>::validParams();
      24       12980 :   MooseEnum reduction("integral average min max");
      25       12980 :   params.addRequiredParam<MooseEnum>("reduction", reduction, "Reduction operation to apply");
      26       12980 :   params.addRequiredParam<std::vector<MooseFunctorName>>("functors",
      27             :                                                          "Functors to apply the reduction on");
      28       12980 :   params.addRequiredParam<MeshDivisionName>(
      29             :       "mesh_division",
      30             :       "Mesh division object which dictates the elements to perform the reduction with");
      31        3245 :   params.addClassDescription("Perform reductions on functors based on a per-mesh-division basis");
      32        6490 :   return params;
      33        3245 : }
      34             : 
      35         100 : MeshDivisionFunctorReductionVectorPostprocessor::MeshDivisionFunctorReductionVectorPostprocessor(
      36         100 :     const InputParameters & parameters)
      37             :   : SpatialUserObjectFunctor<ElementVectorPostprocessor>(parameters),
      38         100 :     _reduction(getParam<MooseEnum>("reduction")),
      39         200 :     _nfunctors(getParam<std::vector<MooseFunctorName>>("functors").size()),
      40         400 :     _mesh_division(_fe_problem.getMeshDivision(getParam<MeshDivisionName>("mesh_division"), _tid))
      41             : {
      42             :   // Gather the functors
      43         200 :   const auto & functor_names = getParam<std::vector<MooseFunctorName>>("functors");
      44         425 :   for (const auto & functor_name : functor_names)
      45         325 :     _functors.push_back(&getFunctor<Real>(functor_name));
      46             : 
      47             :   // Set up reduction vectors
      48         100 :   const auto ndiv = _mesh_division.getNumDivisions();
      49         100 :   _volumes.resize(ndiv);
      50         425 :   for (const auto & functor_name : functor_names)
      51             :   {
      52         325 :     auto & p = declareVector(functor_name);
      53         325 :     p.resize(ndiv);
      54         325 :     _functor_reductions.push_back(&p);
      55             :   }
      56         100 : }
      57             : 
      58             : void
      59          96 : MeshDivisionFunctorReductionVectorPostprocessor::initialize()
      60             : {
      61         408 :   for (auto & reduction : _functor_reductions)
      62             :   {
      63         312 :     if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
      64         192 :       std::fill(reduction->begin(), reduction->end(), 0);
      65         120 :     else if (_reduction == ReductionEnum::MIN)
      66          96 :       std::fill(reduction->begin(), reduction->end(), std::numeric_limits<Real>::max());
      67          24 :     else if (_reduction == ReductionEnum::MAX)
      68          24 :       std::fill(reduction->begin(), reduction->end(), std::numeric_limits<Real>::min());
      69             :     else
      70             :       mooseAssert(false, "Unknown reduction type");
      71             :   }
      72          96 :   std::fill(_volumes.begin(), _volumes.end(), 0);
      73          96 : }
      74             : 
      75             : void
      76       11204 : MeshDivisionFunctorReductionVectorPostprocessor::execute()
      77             : {
      78       11204 :   const auto state_arg = determineState();
      79       11204 :   if (hasBlocks(_current_elem->subdomain_id()))
      80             :   {
      81       11204 :     const auto index = _mesh_division.divisionIndex(*_current_elem);
      82       11204 :     if (index == MooseMeshDivision::INVALID_DIVISION_INDEX)
      83             :     {
      84           8 :       mooseWarning("Spatial value sampled outside of the mesh_division specified in element: " +
      85           8 :                    Moose::stringify(_current_elem->id()) + " of centroid " +
      86           4 :                    Moose::stringify(_current_elem->true_centroid()));
      87           0 :       return;
      88             :     }
      89       47600 :     for (const auto i : make_range(_nfunctors))
      90       36400 :       if (_functors[i]->hasBlocks(_current_elem->subdomain_id()))
      91       72800 :         for (const auto qp : make_range(_qrule->n_points()))
      92             :         {
      93       36400 :           const Moose::ElemQpArg elem_qp = {_current_elem, qp, _qrule, _q_point[qp]};
      94       36400 :           const auto functor_value = (*_functors[i])(elem_qp, state_arg);
      95       36400 :           if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
      96       22400 :             (*_functor_reductions[i])[index] += _JxW[qp] * _coord[qp] * functor_value;
      97       14000 :           else if (_reduction == ReductionEnum::MIN)
      98             :           {
      99       11200 :             if ((*_functor_reductions[i])[index] > functor_value)
     100         660 :               (*_functor_reductions[i])[index] = functor_value;
     101             :           }
     102        2800 :           else if (_reduction == ReductionEnum::MAX)
     103        2800 :             if ((*_functor_reductions[i])[index] < functor_value)
     104        1270 :               (*_functor_reductions[i])[index] = functor_value;
     105             : 
     106       36400 :           if (i == 0 && _reduction == ReductionEnum::AVERAGE)
     107        2800 :             _volumes[index] += _JxW[qp] * _coord[qp];
     108             :         }
     109             :   }
     110             : }
     111             : 
     112             : void
     113          68 : MeshDivisionFunctorReductionVectorPostprocessor::finalize()
     114             : {
     115         289 :   for (auto & reduction : _functor_reductions)
     116         221 :     if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
     117         136 :       gatherSum(*reduction);
     118          85 :     else if (_reduction == ReductionEnum::MIN)
     119          68 :       gatherMin(*reduction);
     120          17 :     else if (_reduction == ReductionEnum::MAX)
     121          17 :       gatherMax(*reduction);
     122             : 
     123          68 :   if (_reduction == ReductionEnum::AVERAGE)
     124             :   {
     125          17 :     gatherSum(_volumes);
     126          85 :     for (const auto i_f : make_range(_nfunctors))
     127        1156 :       for (const auto i : index_range(*_functor_reductions[i_f]))
     128        1088 :         if (!MooseUtils::absoluteFuzzyEqual(_volumes[i], 0))
     129         544 :           (*_functor_reductions[i_f])[i] /= _volumes[i];
     130             :         else
     131         544 :           (*_functor_reductions[i_f])[i] = 0;
     132             :   }
     133          68 : }
     134             : 
     135             : void
     136          12 : MeshDivisionFunctorReductionVectorPostprocessor::threadJoin(const UserObject & s)
     137             : {
     138          12 :   const auto & sibling = static_cast<const MeshDivisionFunctorReductionVectorPostprocessor &>(s);
     139             : 
     140          51 :   for (const auto i_f : make_range(_nfunctors))
     141         663 :     for (const auto i : index_range(*_functor_reductions[i_f]))
     142             :     {
     143         624 :       if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
     144         384 :         (*_functor_reductions[i_f])[i] += (*sibling._functor_reductions[i_f])[i];
     145         240 :       else if (_reduction == ReductionEnum::MIN)
     146             :       {
     147         192 :         if ((*_functor_reductions[i_f])[i] > (*sibling._functor_reductions[i_f])[i])
     148          48 :           (*_functor_reductions[i_f])[i] = (*sibling._functor_reductions[i_f])[i];
     149             :       }
     150          48 :       else if (_reduction == ReductionEnum::MAX)
     151          48 :         if ((*_functor_reductions[i_f])[i] < (*sibling._functor_reductions[i_f])[i])
     152          12 :           (*_functor_reductions[i_f])[i] = (*sibling._functor_reductions[i_f])[i];
     153             : 
     154             :       // Average-reduction requires the reduction of the volume
     155         624 :       if (i_f == 0 && _reduction == ReductionEnum::AVERAGE)
     156          48 :         _volumes[i] += sibling._volumes[i];
     157             :     }
     158          12 : }
     159             : 
     160             : Real
     161          28 : MeshDivisionFunctorReductionVectorPostprocessor::spatialValue(const Point & p) const
     162             : {
     163          28 :   if (_nfunctors > 1)
     164           3 :     mooseError("The spatialValue user object interface was not conceived for objects that compute "
     165             :                "multiple values for a given spatial point. Please specify only one functor");
     166          25 :   const auto index = _mesh_division.divisionIndex(p);
     167          25 :   if (index == MooseMeshDivision::INVALID_DIVISION_INDEX)
     168           3 :     mooseError("Spatial value sampled outside of the mesh_division specified at", p);
     169          22 :   return (*_functor_reductions[0])[index];
     170             : }
     171             : 
     172             : bool
     173       11204 : MeshDivisionFunctorReductionVectorPostprocessor::hasBlocks(const SubdomainID sub) const
     174             : {
     175       11204 :   return BlockRestrictable::hasBlocks(sub);
     176             : }