NekSideIntegralVariableUserObject.C 5.12 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
/****************************************************************/
/*               DO NOT MODIFY THIS HEADER                      */
/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
/*                                                              */
/*           (c) 2010 Battelle Energy Alliance, LLC             */
/*                   ALL RIGHTS RESERVED                        */
/*                                                              */
/*          Prepared by Battelle Energy Alliance, LLC           */
/*            Under Contract No. DE-AC07-05ID14517              */
/*            With the U. S. Department of Energy               */
/*                                                              */
/*            See COPYRIGHT for full restrictions               */
/****************************************************************/

#include "NekSideIntegralVariableUserObject.h"
16 17
#include "MooseVariableScalar.h"
#include "SystemBase.h"
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

template<>
InputParameters validParams<NekSideIntegralVariableUserObject>()
{
  InputParameters params = validParams<SideIntegralUserObject>();
  params.addRequiredCoupledVar("variable", "The name of the variable that will be integrated");
  params.addRequiredParam<std::string>("legendre_function_name", "Name of function to compute Legendre polynomial value at a point");
  params.addRequiredParam<std::string>("fourier_function_name", "Name of function to compute Fourier polynomial value at a point");
  params.addRequiredParam<int>("l_direction", "Direction of integration for Legendre polynomial");
  params.addRequiredParam<int>("l_order", "The order of the Legendre expansion");
  params.addRequiredParam<int>("f_order", "The order of the Fourier expansion");
  params.addRequiredParam<std::string>("aux_scalar_name", "Aux scalar to store the Legendre expansion coefficient");
  params.addParam<MaterialPropertyName>("diffusion_coefficient_name",
                                        "thermal_conductivity",
                                        "Property name of the diffusivity (Default: thermal_conductivity)");
33
  params.addRequiredParam<std::string>("surface_area_pp", "The name of the post processor that calculates surface area");
34 35 36 37 38 39 40
  return params;
}

NekSideIntegralVariableUserObject::NekSideIntegralVariableUserObject(const InputParameters & parameters) :
    // TODO we really shouldn't have to dynamic cast into FourierPolynomial and Legendre Polynomial here
    // but this was a quick example
    SideIntegralUserObject(parameters),
41
    MooseVariableInterface<Real>(this, false),
42 43 44 45 46 47 48 49
    _u(coupledValue("variable")),
    _legendre_poly_function(dynamic_cast<LegendrePolynomial&>(_mci_feproblem.getFunction(parameters.get<std::string>("legendre_function_name")))),
    _fourier_poly_function(dynamic_cast<FourierPolynomial&>(_mci_feproblem.getFunction(parameters.get<std::string>("fourier_function_name")))),
    _l_direction(parameters.get<int>("l_direction")),
    _l_order(parameters.get<int>("l_order")),
    _f_order(parameters.get<int>("f_order")),
    _aux_scalar_name(parameters.get<std::string>("aux_scalar_name")),
    _coupled_grad(coupledGradient("variable")),
50 51
    _diffusion_coefficient(getMaterialProperty<Real>("diffusion_coefficient_name")),
    _surface_area_pp(getPostprocessorValueByName(parameters.get<std::string>("surface_area_pp")))
52

53 54 55 56 57 58 59 60 61 62 63 64 65 66
{
  addMooseVariableDependency(mooseVariable());

  // For now hard code that l_direction better be 2
  // TODO come back and put logic in to determine f diections given l directions
  if ( _l_direction == 2)
  {
    _f_direction_1 = 0;
    _f_direction_2 = 1;
  }
  else
  {
    mooseError("Need to implement logic for l direction not equal to 2");
  }
67 68


69 70 71 72 73 74 75
}

Real
NekSideIntegralVariableUserObject::computeQpIntegral()
{
  Real l_poly_val = _legendre_poly_function.getPolynomialValue(_t,_q_point[_qp](_l_direction),_l_order);
  Real f_poly_val = _fourier_poly_function.getPolynomialValue(_t,_q_point[_qp](_f_direction_1),_q_point[_qp](_f_direction_2), _f_order);
76 77 78
  // There is an added correction factor that accounts for the surface area of the pin
  // The analytic expression is (2 / R / delta z) but we compute numerically as
  // (2 / (Surf Area / (2*PI) ) ) = 4PI / Surf Area
79 80 81 82 83 84 85 86

  //  Implementation from April Novak, to be tested more carefully
  //  RealVectorValue coupled_grad_vec(_coupled_grad[_qp](0), _coupled_grad[_qp](1),  _coupled_grad[_qp](2));
  // RealVectorValue normals_vec(_normals[_qp](0), _normals[_qp](1), _normals[_qp](2));
  // return -_diffusion_coefficient[_qp] * coupled_grad_vec * normals_vec * l_poly_val * f_poly_val * 4.0*M_PI / _surface_area_pp;

  // Original implementation from Matt Ellis
  return -_diffusion_coefficient[_qp] * _coupled_grad[_qp].norm() * l_poly_val * f_poly_val * 4.0*M_PI / _surface_area_pp;
87 88 89 90 91 92 93 94 95 96 97 98
}
void
NekSideIntegralVariableUserObject::finalize()
{
  // In the finalize step store the result of the integral in the scalar variable
  MooseVariableScalar & scalar =  _fe_problem.getScalarVariable(_tid, _aux_scalar_name);
  scalar.reinit();
  // The dof indices for the scalar variable
  std::vector<dof_id_type> & dof = scalar.dofIndices();
  scalar.sys().solution().set(dof[_l_order], getValue());
  scalar.sys().solution().close();
}