Skip to content

Commit

Permalink
adding reporting id vpp, refs idaholab#19217
Browse files Browse the repository at this point in the history
  • Loading branch information
@yjung-anl
yjung-anl committed Nov 6, 2021
1 parent ae7a46d commit 752e5bf
Show file tree
Hide file tree
Showing 5 changed files with 397 additions and 0 deletions.
48 changes: 48 additions & 0 deletions modules/reactor/include/vectorpostprocessors/IntegralReportingIDVPP.h
View file
@@ -0,0 +1,48 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
/*************************************************/

#pragma once

#include "ElementVariableVectorPostprocessor.h"

#include <ostream>

/**
* This IntegralReportingIDVPP source code is to integrate variables based on reporting IDs
*/
class IntegralReportingIDVPP : public ElementVariableVectorPostprocessor
{
public:
static InputParameters validParams();
IntegralReportingIDVPP(const InputParameters & parameters);

virtual void initialize() override;
virtual void execute() override;
virtual void finalize() override;
virtual void threadJoin(const UserObject & uo) override;

protected:
/// Number of variables to be integrated on pins
const unsigned int _nvar;
/// Reporting IDs in use
const std::vector<ExtraElementIDName> _reporting_id;
/// Number of reporting IDs in use
const unsigned int _n_reporting_id;
// Map of element ids to parsed vpp ids
std::map<dof_id_type, dof_id_type> _unique_vpp_ids;
/// Vectors holding reporting IDs
std::vector<VectorPostprocessorValue *> _var_reporting_ids;
/// Coupled MOOSE variables to be integrated
std::vector<const MooseVariable *> _vars;
/// Quadrature point values of coupled MOOSE variables
std::vector<const VariableValue *> _var_values;
/// Vectors holding variable integrals over reporting IDs
std::vector<VectorPostprocessorValue *> _var_integrals;
};
165 changes: 165 additions & 0 deletions modules/reactor/src/vectorpostprocessors/IntegralReportingIDVPP.C
View file
@@ -0,0 +1,165 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "IntegralReportingIDVPP.h"

#include "MooseVariable.h"
#include "MooseMesh.h"

#include "libmesh/mesh_base.h"

registerMooseObject("ReactorApp", IntegralReportingIDVPP);

InputParameters
IntegralReportingIDVPP::validParams()
{
InputParameters params = ElementVariableVectorPostprocessor::validParams();
params.addRequiredParam<std::vector<ExtraElementIDName>>(
"reporting_id", "The types of reporting id by which to separate integrals.");
params.addClassDescription(
"This IntegralReportingIDVPP source code is to integrate variables based on reporting IDs");
return params;
}

IntegralReportingIDVPP::IntegralReportingIDVPP(const InputParameters & parameters)
: ElementVariableVectorPostprocessor(parameters),
_nvar(coupledComponents("variable")),
_reporting_id(getParam<std::vector<ExtraElementIDName>>("reporting_id")),
_n_reporting_id(_reporting_id.size())
{
// numbers of parsed vpp ids
dof_id_type n_vpp_ids;
// list of extra_integer id
std::set<dof_id_type> ids;

std::set<ExtraElementIDName> _reporting_id_copy(_reporting_id.begin(), _reporting_id.end());
// initial setting for parsing
ExtraElementIDName id_type = *_reporting_id_copy.begin();
_reporting_id_copy.erase(id_type);
unsigned int id_index = getElementIDIndexByName(id_type);
if (blockRestricted())
ids = getElemIDsOnBlocks(id_index, blockIDs());
else
ids = getAllElemIDs(id_index);
n_vpp_ids = ids.size();
for (auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
// initializing vpp_id based on the first reporting id value
auto id = elem->get_extra_integer(id_index);
_unique_vpp_ids[elem->id()] = std::distance(ids.begin(), std::find(ids.begin(), ids.end(), id));
}

// looping over reporting ids
while (!_reporting_id_copy.empty())
{
id_type = *_reporting_id_copy.begin();
_reporting_id_copy.erase(id_type);
id_index = getElementIDIndexByName(id_type);
// parsing reporting ids (find unique combinations of reporting ids)
std::set<std::pair<dof_id_type, dof_id_type>> unique_ids;
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
auto elem_id = elem->id();
dof_id_type id1 = _unique_vpp_ids[elem_id];
dof_id_type id2 = elem->get_extra_integer(id_index);
auto it = unique_ids.find(std::pair<dof_id_type, dof_id_type>(id1, id2));
if (it == unique_ids.end())
unique_ids.insert(std::pair<dof_id_type, dof_id_type>(id1, id2));
}
comm().set_union(unique_ids);
// update parsed vpp ids
n_vpp_ids = unique_ids.size();
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
dof_id_type id1 = _unique_vpp_ids[elem->id()];
dof_id_type id2 = elem->get_extra_integer(id_index);
_unique_vpp_ids[elem->id()] = std::distance(
unique_ids.begin(),
std::find(
unique_ids.begin(), unique_ids.end(), std::pair<dof_id_type, dof_id_type>(id1, id2)));
}
}

// set up variable vector containing parsed reporting ids
for (unsigned int i = 0; i < _n_reporting_id; ++i)
{
id_type = _reporting_id[i];
auto & p = declareVector("Level-" + std::to_string(i) + "-" + id_type);
p.resize(n_vpp_ids);
// collect local map of local vpp id to reporting id
std::map<dof_id_type, dof_id_type> reporting_ids;
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
auto vpp_id = _unique_vpp_ids[elem->id()];
if (reporting_ids.find(vpp_id) == reporting_ids.end())
reporting_ids[vpp_id] = elem->get_extra_integer(getElementIDIndexByName(id_type));
}
// create global map of local vpp id to reporting id
comm().set_union(reporting_ids);
for (auto it : reporting_ids)
p[it.first] = it.second;
_var_reporting_ids.push_back(&p);
}

// declare vectors containg integral values
for (unsigned int i = 0; i < _nvar; ++i)
{
_vars.push_back(getVar("variable", i));
_var_values.push_back(&coupledValue("variable", i));
auto & p = declareVector(_vars[i]->name());
p.resize(n_vpp_ids);
_var_integrals.push_back(&p);
}
}

void
IntegralReportingIDVPP::initialize()
{
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
for (size_t i = 0; i < (*_var_integrals[ivar]).size(); ++i)
(*_var_integrals[ivar])[i] = 0;
}

void
IntegralReportingIDVPP::execute()
{
if (hasBlocks(_current_elem->subdomain_id()))
{
auto ipos = _unique_vpp_ids[_current_elem->id()];
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
if (_vars[ivar]->hasBlocks(_current_elem->subdomain_id()))
for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
(*_var_integrals[ivar])[ipos] += _JxW[qp] * _coord[qp] * (*_var_values[ivar])[qp];
}
}

void
IntegralReportingIDVPP::finalize()
{
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
gatherSum(*_var_integrals[ivar]);
}

void
IntegralReportingIDVPP::threadJoin(const UserObject & s)
{
const IntegralReportingIDVPP & sibling = static_cast<const IntegralReportingIDVPP &>(s);

for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
for (size_t i = 0; i < (*_var_integrals[ivar]).size(); ++i)
(*_var_integrals[ivar])[i] += (*sibling._var_integrals[ivar])[i];
}
65 changes: 65 additions & 0 deletions modules/reactor/test/tests/vectorpostprocessors/gold/reporting_id_vpp_out_integral_0001.csv
View file
@@ -0,0 +1,65 @@
Level-0-assembly_id,Level-1-pin_id,value
0,0,1.5671877195416
0,1,1.5069863622713
0,2,1.3888237414075
0,3,1.2173342143126
0,4,1.5069863622713
0,5,1.4490498864433
0,6,1.3354739906409
0,7,1.1705334893373
0,8,1.3888237414075
0,9,1.3354739906409
0,10,1.2307596343731
0,11,1.0787875868243
0,12,1.2173342143126
0,13,1.1705334893373
0,14,1.0787875868243
0,15,0.9455496753278
1,0,0.99904037162165
1,1,0.74234175714638
1,2,0.45713179810106
1,3,0.15435705940956
1,4,0.96063200921654
1,5,0.71382572120417
1,6,0.43957170982704
1,7,0.14842276281255
1,8,0.88532353029201
1,9,0.657865552518
1,10,0.40511160801102
1,11,0.13678720164245
1,12,0.77600532884632
1,13,0.5766144950035
1,14,0.35507745371261
1,15,0.11989695716944
2,0,0.99904037162165
2,1,0.96063200921654
2,2,0.88532353029201
2,3,0.77600532884632
2,4,0.74234175714638
2,5,0.71382572120417
2,6,0.657865552518
2,7,0.5766144950035
2,8,0.45713179810106
2,9,0.43957170982704
2,10,0.40511160801102
2,11,0.35507745371261
2,12,0.15435705940956
2,13,0.14842276281255
2,14,0.13678720164245
2,15,0.11989695716944
3,0,0.63684064003929
3,1,0.47322307065213
3,2,0.29140453553409
3,3,0.098396889891493
3,4,0.47322307065213
3,5,0.35163068056988
3,6,0.21653666621981
3,7,0.073114282680312
3,8,0.29140453553409
3,9,0.21653666621981
3,10,0.1333404277161
3,11,0.045024293667973
3,12,0.098396889891493
3,13,0.073114282680312
3,14,0.045024293667973
3,15,0.015202593593918
108 changes: 108 additions & 0 deletions modules/reactor/test/tests/vectorpostprocessors/reporting_id_vpp.i
View file
@@ -0,0 +1,108 @@
[Mesh]
[pin1]
type = ConcentricCircleMeshGenerator
num_sectors = 4
radii = '0.4 0.5'
rings = '2 1 1'
has_outer_square = on
pitch = 1.26
preserve_volumes = yes
smoothing_max_it = 3
[]

[pin2]
type = ConcentricCircleMeshGenerator
num_sectors = 4
radii = '0.3 0.4'
rings = '2 1 1'
has_outer_square = on
pitch = 1.26
preserve_volumes = yes
smoothing_max_it = 3
[]

[assembly1]
type = CartesianIDPatternedMeshGenerator
inputs = 'pin1 pin2'
pattern = ' 1 0 1 0;
0 1 0 1;
1 0 1 0;
0 1 0 1'
assign_type = 'cell'
id_name = 'pin_id'
[]

[assembly2]
type = CartesianIDPatternedMeshGenerator
inputs = 'pin1 pin2'
pattern = ' 0 1 1 0;
1 0 0 1;
1 0 0 1;
0 1 1 0'
assign_type = 'cell'
id_name = 'pin_id'
[]

[core]
type = CartesianIDPatternedMeshGenerator
inputs = 'assembly1 assembly2'
pattern = '0 1;
1 0'
assign_type = 'cell'
id_name = 'assembly_id'
[]
[]

[Executioner]
type = Steady
[]

[Problem]
solve = false
[]

[AuxVariables]
[pin_id]
family = MONOMIAL
order = CONSTANT
[]
[assembly_id]
family = MONOMIAL
order = CONSTANT
[]
[value]
order = FIRST
[]
[]

[AuxKernels]
[set_pin_id]
type = ElemExtraIDAux
variable = pin_id
extra_id_name = pin_id
[]
[set_assembly_id]
type = ElemExtraIDAux
variable = assembly_id
extra_id_name = assembly_id
[]
[set_value]
type = ParsedAux
use_xyzt = true
function = 'cos(3.14159265*(x+0.63)/(20.16)) * cos(3.14159265*(-y+0.63)/(20.16))'
variable = value
[]
[]

[VectorPostprocessors]
[integral]
type = IntegralReportingIDVPP
variable = 'value'
reporting_id = 'assembly_id pin_id'
[]
[]
[Outputs]
exodus = false
csv = true
execute_on = timestep_end
[]

0 comments on commit 752e5bf

Please sign in to comment.