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 "ExtraIDIntegralVectorPostprocessor.h"
      11             : #include "MooseVariable.h"
      12             : #include "MooseMesh.h"
      13             : #include "MooseMeshUtils.h"
      14             : 
      15             : #include "libmesh/mesh_base.h"
      16             : 
      17             : registerMooseObject("MooseApp", ExtraIDIntegralVectorPostprocessor);
      18             : 
      19             : InputParameters
      20       14494 : ExtraIDIntegralVectorPostprocessor::validParams()
      21             : {
      22       14494 :   InputParameters params = ElementVariableVectorPostprocessor::validParams();
      23       57976 :   params.addParam<std::vector<MaterialPropertyName>>(
      24             :       "mat_prop", "The names of material properties that this VectorPostprocessor operates on");
      25       57976 :   params.addRequiredParam<std::vector<ExtraElementIDName>>(
      26             :       "id_name", "List of extra element ID names by which to separate integral(s).");
      27       57976 :   params.addParam<bool>("average", false, "Whether or not to compute volume average");
      28       57976 :   params.addParam<std::string>("spatial_value_name",
      29             :                                "To specify which variable or material property is to be used when "
      30             :                                "using as a spatial user object functor");
      31       28988 :   params.addClassDescription("Integrates or averages variables based on extra element IDs");
      32       14494 :   params.makeParamNotRequired("variable");
      33       14494 :   return params;
      34           0 : }
      35             : 
      36         120 : ExtraIDIntegralVectorPostprocessor::ExtraIDIntegralVectorPostprocessor(
      37         120 :     const InputParameters & parameters)
      38             :   : SpatialUserObjectFunctor<ElementVariableVectorPostprocessor>(parameters),
      39         120 :     _average(getParam<bool>("average")),
      40         480 :     _nvar(isParamValid("variable") ? coupledComponents("variable") : 0),
      41         288 :     _nprop(isParamValid("mat_prop") ? getParam<std::vector<MaterialPropertyName>>("mat_prop").size()
      42             :                                     : 0),
      43         288 :     _prop_names(isParamValid("mat_prop") ? getParam<std::vector<MaterialPropertyName>>("mat_prop")
      44             :                                          : std::vector<MaterialPropertyName>()),
      45         240 :     _extra_id(getParam<std::vector<ExtraElementIDName>>("id_name")),
      46         120 :     _n_extra_id(_extra_id.size()),
      47         240 :     _spatial_evaluation_index(std::numeric_limits<unsigned int>::max())
      48             : {
      49         120 :   if (!_nvar && !_nprop)
      50           0 :     mooseError("Neither 'variable' nor 'mat_prop' was specified.");
      51             : 
      52             :   // create map of element ids to parsed vpp ids
      53             :   _unique_vpp_ids =
      54         120 :       MooseMeshUtils::getExtraIDUniqueCombinationMap(_mesh.getMesh(), blockIDs(), _extra_id);
      55             : 
      56             :   // set up variable vector containing parsed extra ids
      57         312 :   for (unsigned int i = 0; i < _n_extra_id; ++i)
      58             :   {
      59         192 :     const auto id_type = _extra_id[i];
      60         192 :     auto & p = declareVector("Level-" + std::to_string(i) + "-" + id_type);
      61             :     // collect local map of local vpp id to extra id
      62         192 :     std::map<dof_id_type, dof_id_type> extra_ids;
      63      774861 :     for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
      64             :     {
      65      774669 :       if (!hasBlocks(elem->subdomain_id()))
      66          50 :         continue;
      67      774619 :       auto vpp_id = _unique_vpp_ids[elem->id()];
      68      774619 :       if (extra_ids.find(vpp_id) == extra_ids.end())
      69        7233 :         extra_ids[vpp_id] = elem->get_extra_integer(getElementIDIndexByName(id_type));
      70         192 :     }
      71             :     // create global map of local vpp id to extra id
      72         192 :     comm().set_union(extra_ids);
      73         192 :     p.resize(extra_ids.size());
      74        9552 :     for (auto it : extra_ids)
      75        9360 :       p[it.first] = it.second;
      76         192 :     _extra_ids.push_back(&p);
      77         192 :   }
      78             : 
      79         120 :   std::vector<std::string> vector_names;
      80             : 
      81             :   // declare vectors containing variable integral values
      82         288 :   for (unsigned int i = 0; i < _nvar; ++i)
      83             :   {
      84         336 :     _vars.push_back(getVar("variable", i));
      85         336 :     _var_values.push_back(&coupledValue("variable", i));
      86         168 :     vector_names.push_back(_vars[i]->name());
      87         168 :     auto & p = declareVector(vector_names.back());
      88         168 :     p.resize((*_extra_ids[0]).size());
      89         168 :     _integrals.push_back(&p);
      90             :   }
      91             : 
      92             :   // declare vectors containing material property integral values
      93         168 :   for (auto & name : _prop_names)
      94             :   {
      95          48 :     _props.push_back(&getMaterialPropertyByName<Real>(name));
      96          48 :     vector_names.push_back(name);
      97          48 :     auto & p = declareVector(vector_names.back());
      98          48 :     p.resize((*_extra_ids[0]).size());
      99          48 :     _integrals.push_back(&p);
     100             :   }
     101             : 
     102         360 :   if (isParamValid("spatial_value_name"))
     103             :   {
     104          38 :     const auto & name = getParam<std::string>("spatial_value_name");
     105          19 :     _spatial_evaluation_index = vector_names.size();
     106          19 :     for (const auto i : index_range(vector_names))
     107          19 :       if (name == vector_names[i])
     108             :       {
     109          19 :         _spatial_evaluation_index = i;
     110          19 :         break;
     111             :       }
     112          19 :     if (_spatial_evaluation_index == vector_names.size())
     113           0 :       paramError("spatial_value_name",
     114             :                  "not a variable name or material property name of this vector postprocessor");
     115             :   }
     116         120 : }
     117             : 
     118             : void
     119         119 : ExtraIDIntegralVectorPostprocessor::initialize()
     120             : {
     121         334 :   for (auto & integral : _integrals)
     122         215 :     std::fill(integral->begin(), integral->end(), 0);
     123             : 
     124         119 :   if (_average)
     125          12 :     _volumes.assign((*_extra_ids[0]).size(), 0);
     126         119 : }
     127             : 
     128             : void
     129      393601 : ExtraIDIntegralVectorPostprocessor::execute()
     130             : {
     131      393601 :   if (hasBlocks(_current_elem->subdomain_id()))
     132             :   {
     133      393601 :     unsigned int i = 0;
     134      393601 :     auto ipos = _unique_vpp_ids[_current_elem->id()];
     135      983810 :     for (unsigned int ivar = 0; ivar < _nvar; ++ivar, ++i)
     136      590209 :       if (_vars[ivar]->hasBlocks(_current_elem->subdomain_id()))
     137     1181573 :         for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     138      591364 :           (*_integrals[i])[ipos] += _JxW[qp] * _coord[qp] * (*_var_values[ivar])[qp];
     139             : 
     140      590209 :     for (unsigned int iprop = 0; iprop < _nprop; ++iprop, ++i)
     141      393216 :       for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     142      196608 :         (*_integrals[i])[ipos] += _JxW[qp] * _coord[qp] * (*_props[iprop])[qp];
     143             : 
     144      393601 :     if (_average)
     145       49152 :       _volumes[ipos] += _current_elem->volume();
     146             :   }
     147      393601 : }
     148             : 
     149             : void
     150         108 : ExtraIDIntegralVectorPostprocessor::finalize()
     151             : {
     152         304 :   for (auto & integral : _integrals)
     153         196 :     gatherSum(*integral);
     154             : 
     155         108 :   if (_average)
     156             :   {
     157          11 :     gatherSum(_volumes);
     158             : 
     159          55 :     for (auto & integral : _integrals)
     160        2860 :       for (unsigned int i = 0; i < integral->size(); ++i)
     161        2816 :         (*integral)[i] /= _volumes[i];
     162             :   }
     163         108 : }
     164             : 
     165             : void
     166          11 : ExtraIDIntegralVectorPostprocessor::threadJoin(const UserObject & s)
     167             : {
     168          11 :   const auto & sibling = static_cast<const ExtraIDIntegralVectorPostprocessor &>(s);
     169             : 
     170          30 :   for (unsigned int i = 0; i < _integrals.size(); ++i)
     171         929 :     for (size_t j = 0; j < (*_integrals[i]).size(); ++j)
     172         910 :       (*_integrals[i])[j] += (*sibling._integrals[i])[j];
     173             : 
     174          11 :   if (_average)
     175          65 :     for (unsigned int i = 0; i < _volumes.size(); ++i)
     176          64 :       _volumes[i] += sibling._volumes[i];
     177          11 : }
     178             : 
     179             : Real
     180           0 : ExtraIDIntegralVectorPostprocessor::evaluate(const ElemArg & elem,
     181             :                                              const Moose::StateArg & libmesh_dbg_var(state)) const
     182             : {
     183             :   mooseAssert(state.state == 0, "We do not currently support evaluating at old states");
     184           0 :   return elementValue(elem.elem);
     185             : }
     186             : 
     187             : Real
     188         961 : ExtraIDIntegralVectorPostprocessor::evaluate(const ElemQpArg & elem,
     189             :                                              const Moose::StateArg & libmesh_dbg_var(state)) const
     190             : {
     191             :   mooseAssert(state.state == 0, "We do not currently support evaluating at old states");
     192         961 :   return elementValue(elem.elem);
     193             : }
     194             : 
     195             : Real
     196         961 : ExtraIDIntegralVectorPostprocessor::elementValue(const Elem * elem) const
     197             : {
     198         961 :   if (_spatial_evaluation_index == std::numeric_limits<unsigned int>::max())
     199          10 :     paramError("spatial_value_name",
     200             :                "Must set when ExtraIDIntegralVectorPostprocessor is used as a functor");
     201         956 :   const auto it = _unique_vpp_ids.find(elem->id());
     202         956 :   if (it == _unique_vpp_ids.end())
     203             :   {
     204           4 :     if (!hasBlocks(elem->subdomain_id()))
     205           8 :       mooseError(
     206             :           "Failed evaluating spatial value of element ",
     207           4 :           elem->id(),
     208             :           " in subdomain ",
     209           4 :           elem->subdomain_id(),
     210             :           ". This is likely due to the object using this "
     211             :           "ExtraIDIntegralVectorPostprocessor as a functor operates on a subdomain outside of "
     212             :           "this ExtraIDIntegralVectorPostprocessor");
     213             :     else
     214           0 :       mooseError("Failed evaluating spatial value of element ",
     215           0 :                  elem->id(),
     216             :                  ". We should not hit this error. Please contact code developers for help");
     217             :   }
     218        1904 :   return (*_integrals[_spatial_evaluation_index])[it->second];
     219             : }