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       14481 : MeshDivisionFunctorReductionVectorPostprocessor::validParams()
      22             : {
      23       14481 :   InputParameters params = SpatialUserObjectFunctor<ElementVectorPostprocessor>::validParams();
      24       57924 :   MooseEnum reduction("integral average min max");
      25       57924 :   params.addRequiredParam<MooseEnum>("reduction", reduction, "Reduction operation to apply");
      26       57924 :   params.addRequiredParam<std::vector<MooseFunctorName>>("functors",
      27             :                                                          "Functors to apply the reduction on");
      28       57924 :   params.addRequiredParam<MeshDivisionName>(
      29             :       "mesh_division",
      30             :       "Mesh division object which dictates the elements to perform the reduction with");
      31       14481 :   params.addClassDescription("Perform reductions on functors based on a per-mesh-division basis");
      32       28962 :   return params;
      33       14481 : }
      34             : 
      35         116 : MeshDivisionFunctorReductionVectorPostprocessor::MeshDivisionFunctorReductionVectorPostprocessor(
      36         116 :     const InputParameters & parameters)
      37             :   : SpatialUserObjectFunctor<ElementVectorPostprocessor>(parameters),
      38         116 :     _reduction(getParam<MooseEnum>("reduction")),
      39         232 :     _nfunctors(getParam<std::vector<MooseFunctorName>>("functors").size()),
      40         464 :     _mesh_division(_fe_problem.getMeshDivision(getParam<MeshDivisionName>("mesh_division"), _tid))
      41             : {
      42             :   // Gather the functors
      43         232 :   const auto & functor_names = getParam<std::vector<MooseFunctorName>>("functors");
      44         493 :   for (const auto & functor_name : functor_names)
      45         377 :     _functors.push_back(&getFunctor<Real>(functor_name));
      46             : 
      47             :   // Set up reduction vectors
      48         116 :   const auto ndiv = _mesh_division.getNumDivisions();
      49         116 :   _volumes.resize(ndiv);
      50         493 :   for (const auto & functor_name : functor_names)
      51             :   {
      52         377 :     auto & p = declareVector(functor_name);
      53         377 :     p.resize(ndiv);
      54         377 :     _functor_reductions.push_back(&p);
      55             :   }
      56         116 : }
      57             : 
      58             : void
      59         112 : MeshDivisionFunctorReductionVectorPostprocessor::initialize()
      60             : {
      61         476 :   for (auto & reduction : _functor_reductions)
      62             :   {
      63         364 :     if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
      64         224 :       std::fill(reduction->begin(), reduction->end(), 0);
      65         140 :     else if (_reduction == ReductionEnum::MIN)
      66         112 :       std::fill(reduction->begin(), reduction->end(), std::numeric_limits<Real>::max());
      67          28 :     else if (_reduction == ReductionEnum::MAX)
      68          28 :       std::fill(reduction->begin(), reduction->end(), std::numeric_limits<Real>::min());
      69             :     else
      70             :       mooseAssert(false, "Unknown reduction type");
      71             :   }
      72         112 :   std::fill(_volumes.begin(), _volumes.end(), 0);
      73         112 : }
      74             : 
      75             : void
      76       13605 : MeshDivisionFunctorReductionVectorPostprocessor::execute()
      77             : {
      78       13605 :   const auto state_arg = determineState();
      79       13605 :   if (hasBlocks(_current_elem->subdomain_id()))
      80             :   {
      81       13605 :     const auto index = _mesh_division.divisionIndex(*_current_elem);
      82       13605 :     if (index == MooseMeshDivision::INVALID_DIVISION_INDEX)
      83             :     {
      84          10 :       mooseWarning("Spatial value sampled outside of the mesh_division specified in element: " +
      85          10 :                    Moose::stringify(_current_elem->id()) + " of centroid " +
      86           5 :                    Moose::stringify(_current_elem->true_centroid()));
      87           0 :       return;
      88             :     }
      89       57800 :     for (const auto i : make_range(_nfunctors))
      90       44200 :       if (_functors[i]->hasBlocks(_current_elem->subdomain_id()))
      91       88400 :         for (const auto qp : make_range(_qrule->n_points()))
      92             :         {
      93       44200 :           const Moose::ElemQpArg elem_qp = {_current_elem, qp, _qrule, _q_point[qp]};
      94       44200 :           const auto functor_value = (*_functors[i])(elem_qp, state_arg);
      95       44200 :           if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
      96       27200 :             (*_functor_reductions[i])[index] += _JxW[qp] * _coord[qp] * functor_value;
      97       17000 :           else if (_reduction == ReductionEnum::MIN)
      98             :           {
      99       13600 :             if ((*_functor_reductions[i])[index] > functor_value)
     100         798 :               (*_functor_reductions[i])[index] = functor_value;
     101             :           }
     102        3400 :           else if (_reduction == ReductionEnum::MAX)
     103        3400 :             if ((*_functor_reductions[i])[index] < functor_value)
     104        1540 :               (*_functor_reductions[i])[index] = functor_value;
     105             : 
     106       44200 :           if (i == 0 && _reduction == ReductionEnum::AVERAGE)
     107        3400 :             _volumes[index] += _JxW[qp] * _coord[qp];
     108             :         }
     109             :   }
     110             : }
     111             : 
     112             : void
     113          80 : MeshDivisionFunctorReductionVectorPostprocessor::finalize()
     114             : {
     115         340 :   for (auto & reduction : _functor_reductions)
     116         260 :     if (_reduction == ReductionEnum::INTEGRAL || _reduction == ReductionEnum::AVERAGE)
     117         160 :       gatherSum(*reduction);
     118         100 :     else if (_reduction == ReductionEnum::MIN)
     119          80 :       gatherMin(*reduction);
     120          20 :     else if (_reduction == ReductionEnum::MAX)
     121          20 :       gatherMax(*reduction);
     122             : 
     123          80 :   if (_reduction == ReductionEnum::AVERAGE)
     124             :   {
     125          20 :     gatherSum(_volumes);
     126         100 :     for (const auto i_f : make_range(_nfunctors))
     127        1360 :       for (const auto i : index_range(*_functor_reductions[i_f]))
     128        1280 :         if (!MooseUtils::absoluteFuzzyEqual(_volumes[i], 0))
     129         640 :           (*_functor_reductions[i_f])[i] /= _volumes[i];
     130             :         else
     131         640 :           (*_functor_reductions[i_f])[i] = 0;
     132             :   }
     133          80 : }
     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          32 : MeshDivisionFunctorReductionVectorPostprocessor::spatialValue(const Point & p) const
     162             : {
     163          32 :   if (_nfunctors > 1)
     164           4 :     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          28 :   const auto index = _mesh_division.divisionIndex(p);
     167          28 :   if (index == MooseMeshDivision::INVALID_DIVISION_INDEX)
     168           4 :     mooseError("Spatial value sampled outside of the mesh_division specified at", p);
     169          24 :   return (*_functor_reductions[0])[index];
     170             : }
     171             : 
     172             : bool
     173       13605 : MeshDivisionFunctorReductionVectorPostprocessor::hasBlocks(const SubdomainID sub) const
     174             : {
     175       13605 :   return BlockRestrictable::hasBlocks(sub);
     176             : }