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