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       14473 : ExtraIDIntegralVectorPostprocessor::validParams()
      21             : {
      22       14473 :   InputParameters params = ElementVariableVectorPostprocessor::validParams();
      23       14473 :   params.addParam<std::vector<MaterialPropertyName>>(
      24             :       "mat_prop", "The names of material properties that this VectorPostprocessor operates on");
      25       14473 :   params.addRequiredParam<std::vector<ExtraElementIDName>>(
      26             :       "id_name", "List of extra element ID names by which to separate integral(s).");
      27       14473 :   params.addParam<bool>("average", false, "Whether or not to compute volume average");
      28       14473 :   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       14473 :   params.addClassDescription("Integrates or averages variables based on extra element IDs");
      32       14473 :   params.makeParamNotRequired("variable");
      33       14473 :   return params;
      34           0 : }
      35             : 
      36         109 : ExtraIDIntegralVectorPostprocessor::ExtraIDIntegralVectorPostprocessor(
      37         109 :     const InputParameters & parameters)
      38             :   : SpatialUserObjectFunctor<ElementVariableVectorPostprocessor>(parameters),
      39         109 :     _average(getParam<bool>("average")),
      40         109 :     _nvar(isParamValid("variable") ? coupledComponents("variable") : 0),
      41         109 :     _nprop(isParamValid("mat_prop") ? getParam<std::vector<MaterialPropertyName>>("mat_prop").size()
      42             :                                     : 0),
      43         109 :     _prop_names(isParamValid("mat_prop") ? getParam<std::vector<MaterialPropertyName>>("mat_prop")
      44             :                                          : std::vector<MaterialPropertyName>()),
      45         109 :     _extra_id(getParam<std::vector<ExtraElementIDName>>("id_name")),
      46         109 :     _n_extra_id(_extra_id.size()),
      47         218 :     _spatial_evaluation_index(std::numeric_limits<unsigned int>::max())
      48             : {
      49         109 :   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         109 :       MooseMeshUtils::getExtraIDUniqueCombinationMap(_mesh.getMesh(), blockIDs(), _extra_id);
      55             : 
      56             :   // set up variable vector containing parsed extra ids
      57         284 :   for (unsigned int i = 0; i < _n_extra_id; ++i)
      58             :   {
      59         175 :     const auto id_type = _extra_id[i];
      60         175 :     auto & p = declareVector("Level-" + std::to_string(i) + "-" + id_type);
      61             :     // collect local map of local vpp id to extra id
      62         175 :     std::map<dof_id_type, dof_id_type> extra_ids;
      63     1377331 :     for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
      64             :     {
      65      688578 :       if (!hasBlocks(elem->subdomain_id()))
      66          50 :         continue;
      67      688528 :       auto vpp_id = _unique_vpp_ids[elem->id()];
      68      688528 :       if (extra_ids.find(vpp_id) == extra_ids.end())
      69        6451 :         extra_ids[vpp_id] = elem->get_extra_integer(getElementIDIndexByName(id_type));
      70         175 :     }
      71             :     // create global map of local vpp id to extra id
      72         175 :     comm().set_union(extra_ids);
      73         175 :     p.resize(extra_ids.size());
      74        8753 :     for (auto it : extra_ids)
      75        8578 :       p[it.first] = it.second;
      76         175 :     _extra_ids.push_back(&p);
      77         175 :   }
      78             : 
      79         109 :   std::vector<std::string> vector_names;
      80             : 
      81             :   // declare vectors containing variable integral values
      82         262 :   for (unsigned int i = 0; i < _nvar; ++i)
      83             :   {
      84         153 :     _vars.push_back(getVar("variable", i));
      85         153 :     _var_values.push_back(&coupledValue("variable", i));
      86         153 :     vector_names.push_back(_vars[i]->name());
      87         153 :     auto & p = declareVector(vector_names.back());
      88         153 :     p.resize((*_extra_ids[0]).size());
      89         153 :     _integrals.push_back(&p);
      90             :   }
      91             : 
      92             :   // declare vectors containing material property integral values
      93         153 :   for (auto & name : _prop_names)
      94             :   {
      95          44 :     _props.push_back(&getMaterialPropertyByName<Real>(name));
      96          44 :     vector_names.push_back(name);
      97          44 :     auto & p = declareVector(vector_names.back());
      98          44 :     p.resize((*_extra_ids[0]).size());
      99          44 :     _integrals.push_back(&p);
     100             :   }
     101             : 
     102         109 :   if (isParamValid("spatial_value_name"))
     103             :   {
     104          18 :     const auto & name = getParam<std::string>("spatial_value_name");
     105          18 :     _spatial_evaluation_index = vector_names.size();
     106          18 :     for (const auto i : index_range(vector_names))
     107          18 :       if (name == vector_names[i])
     108             :       {
     109          18 :         _spatial_evaluation_index = i;
     110          18 :         break;
     111             :       }
     112          18 :     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         109 : }
     117             : 
     118             : void
     119         108 : ExtraIDIntegralVectorPostprocessor::initialize()
     120             : {
     121         304 :   for (auto & integral : _integrals)
     122         196 :     std::fill(integral->begin(), integral->end(), 0);
     123             : 
     124         108 :   if (_average)
     125          11 :     _volumes.assign((*_extra_ids[0]).size(), 0);
     126         108 : }
     127             : 
     128             : void
     129      344399 : ExtraIDIntegralVectorPostprocessor::execute()
     130             : {
     131      344399 :   if (hasBlocks(_current_elem->subdomain_id()))
     132             :   {
     133      344399 :     unsigned int i = 0;
     134      344399 :     auto ipos = _unique_vpp_ids[_current_elem->id()];
     135      860830 :     for (unsigned int ivar = 0; ivar < _nvar; ++ivar, ++i)
     136      516431 :       if (_vars[ivar]->hasBlocks(_current_elem->subdomain_id()))
     137     1033867 :         for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     138      517436 :           (*_integrals[i])[ipos] += _JxW[qp] * _coord[qp] * (*_var_values[ivar])[qp];
     139             : 
     140      516431 :     for (unsigned int iprop = 0; iprop < _nprop; ++iprop, ++i)
     141      344064 :       for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
     142      172032 :         (*_integrals[i])[ipos] += _JxW[qp] * _coord[qp] * (*_props[iprop])[qp];
     143             : 
     144      344399 :     if (_average)
     145       43008 :       _volumes[ipos] += _current_elem->volume();
     146             :   }
     147      344399 : }
     148             : 
     149             : void
     150          98 : ExtraIDIntegralVectorPostprocessor::finalize()
     151             : {
     152         276 :   for (auto & integral : _integrals)
     153         178 :     gatherSum(*integral);
     154             : 
     155          98 :   if (_average)
     156             :   {
     157          10 :     gatherSum(_volumes);
     158             : 
     159          50 :     for (auto & integral : _integrals)
     160        2600 :       for (unsigned int i = 0; i < integral->size(); ++i)
     161        2560 :         (*integral)[i] /= _volumes[i];
     162             :   }
     163          98 : }
     164             : 
     165             : void
     166          10 : ExtraIDIntegralVectorPostprocessor::threadJoin(const UserObject & s)
     167             : {
     168          10 :   const auto & sibling = static_cast<const ExtraIDIntegralVectorPostprocessor &>(s);
     169             : 
     170          28 :   for (unsigned int i = 0; i < _integrals.size(); ++i)
     171         926 :     for (size_t j = 0; j < (*_integrals[i]).size(); ++j)
     172         908 :       (*_integrals[i])[j] += (*sibling._integrals[i])[j];
     173             : 
     174          10 :   if (_average)
     175          65 :     for (unsigned int i = 0; i < _volumes.size(); ++i)
     176          64 :       _volumes[i] += sibling._volumes[i];
     177          10 : }
     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         859 : 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         859 :   return elementValue(elem.elem);
     193             : }
     194             : 
     195             : Real
     196         859 : ExtraIDIntegralVectorPostprocessor::elementValue(const Elem * elem) const
     197             : {
     198         859 :   if (_spatial_evaluation_index == std::numeric_limits<unsigned int>::max())
     199           3 :     paramError("spatial_value_name",
     200             :                "Must set when ExtraIDIntegralVectorPostprocessor is used as a functor");
     201         856 :   const auto it = _unique_vpp_ids.find(elem->id());
     202         856 :   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         852 :   return (*_integrals[_spatial_evaluation_index])[it->second];
     219             : }