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 "UnobstructedPlanarViewFactor.h"
      11             : #include "Assembly.h"
      12             : #include "MathUtils.h"
      13             : 
      14             : #include "libmesh/quadrature.h"
      15             : #include "libmesh/fe_base.h"
      16             : #include "libmesh/mesh_generation.h"
      17             : #include "libmesh/mesh.h"
      18             : #include "libmesh/string_to_enum.h"
      19             : #include "libmesh/quadrature_gauss.h"
      20             : #include "libmesh/point_locator_base.h"
      21             : #include "libmesh/elem.h"
      22             : 
      23             : registerMooseObject("HeatTransferApp", UnobstructedPlanarViewFactor);
      24             : 
      25             : InputParameters
      26         207 : UnobstructedPlanarViewFactor::validParams()
      27             : {
      28         207 :   InputParameters params = ViewFactorBase::validParams();
      29         207 :   params.addClassDescription(
      30             :       "Computes the view factors for planar faces in unubstructed radiative heat transfer.");
      31         207 :   return params;
      32           0 : }
      33             : 
      34         108 : UnobstructedPlanarViewFactor::UnobstructedPlanarViewFactor(const InputParameters & parameters)
      35             :   : ViewFactorBase(parameters),
      36         108 :     _boundary_info(nullptr),
      37             :     _current_remote_side(nullptr),
      38             :     _current_remote_fe(nullptr),
      39         108 :     _current_remote_JxW(nullptr),
      40         108 :     _current_remote_xyz(nullptr),
      41         108 :     _current_remote_normals(nullptr)
      42             : {
      43         108 :   _mesh.errorIfDistributedMesh("UnobstructedPlanarViewFactor");
      44             : 
      45         108 :   if (_mesh.dimension() == 1)
      46           0 :     mooseError("View factor calculations for 1D geometry makes no sense");
      47         108 :   else if (_mesh.dimension() == 2)
      48             :   {
      49          72 :     _exponent = 1;
      50          72 :     _divisor = 2;
      51             :   }
      52             :   else
      53             :   {
      54          36 :     _exponent = 2;
      55          36 :     _divisor = libMesh::pi;
      56             :   }
      57         108 : }
      58             : 
      59             : void
      60        1752 : UnobstructedPlanarViewFactor::execute()
      61             : {
      62        1752 :   auto current_boundary_name = _mesh.getBoundaryName(_current_boundary_id);
      63        1752 :   if (_side_name_index.find(current_boundary_name) == _side_name_index.end())
      64           0 :     mooseError("Current boundary name: ",
      65             :                current_boundary_name,
      66             :                " with id ",
      67           0 :                _current_boundary_id,
      68             :                " not in boundary parameter.");
      69        1752 :   unsigned int index = _side_name_index.find(current_boundary_name)->second;
      70             : 
      71        1752 :   _areas[index] += _current_side_volume;
      72             : 
      73      127800 :   for (auto & side : _side_list)
      74             :   {
      75      126048 :     auto remote_boundary_name = _mesh.getBoundaryName(std::get<2>(side));
      76      126048 :     if (_side_name_index.find(remote_boundary_name) != _side_name_index.end() &&
      77      107040 :         std::get<2>(side) != _current_boundary_id)
      78             :     {
      79             :       // this is the remote side index
      80       86112 :       unsigned int remote_index = _side_name_index.find(remote_boundary_name)->second;
      81       86112 :       Real & vf = _view_factors[index][remote_index];
      82             : 
      83             :       // compute some important quantities on the face
      84       86112 :       reinitFace(std::get<0>(side), std::get<1>(side));
      85             : 
      86             :       // loop over pairs of the qps on the current side
      87      362016 :       for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
      88             :       {
      89             :         // loop over the qps on the remote element
      90     1242432 :         for (unsigned int r_qp = 0; r_qp < _current_remote_JxW->size(); ++r_qp)
      91             :         {
      92      966528 :           Point r2r = (_q_point[qp] - (*_current_remote_xyz)[r_qp]);
      93      966528 :           Real distance = r2r.norm();
      94      966528 :           Real cos1 = r2r * _normals[qp] / distance;
      95      966528 :           Real cos2 = r2r * (*_current_remote_normals)[r_qp] / distance;
      96             : 
      97      966528 :           vf += _JxW[qp] * _coord[qp] * (*_current_remote_JxW)[r_qp] * _current_remote_coord[r_qp] *
      98      966528 :                 std::abs(cos1) * std::abs(cos2) / MathUtils::pow(distance, _exponent);
      99             :         }
     100             :       }
     101             :     }
     102             :   }
     103        1752 : }
     104             : 
     105             : void
     106          81 : UnobstructedPlanarViewFactor::initialize()
     107             : {
     108             :   // get boundary info from the mesh
     109          81 :   _boundary_info = &_mesh.getMesh().get_boundary_info();
     110             : 
     111             :   // get a list of all sides
     112          81 :   _side_list = _boundary_info->build_active_side_list();
     113             : 
     114             :   // set view_factors to zero
     115         639 :   for (unsigned int j = 0; j < _n_sides; ++j)
     116             :   {
     117         558 :     _areas[j] = 0;
     118        4734 :     for (auto & vf : _view_factors[j])
     119        4176 :       vf = 0;
     120             :   }
     121          81 : }
     122             : 
     123             : void
     124          72 : UnobstructedPlanarViewFactor::finalizeViewFactor()
     125             : {
     126          72 :   gatherSum(_areas);
     127             : 
     128             :   // divide view_factor Fij by Ai and pi
     129         568 :   for (unsigned int i = 0; i < _n_sides; ++i)
     130        4208 :     for (auto & vf : _view_factors[i])
     131        3712 :       vf /= (_areas[i] * _divisor);
     132          72 : }
     133             : 
     134             : void
     135           9 : UnobstructedPlanarViewFactor::threadJoinViewFactor(const UserObject & y)
     136             : {
     137             :   const auto & pps = static_cast<const UnobstructedPlanarViewFactor &>(y);
     138          71 :   for (unsigned int i = 0; i < _n_sides; ++i)
     139          62 :     _areas[i] += pps._areas[i];
     140           9 : }
     141             : 
     142             : void
     143       86112 : UnobstructedPlanarViewFactor::reinitFace(dof_id_type elem_id, unsigned int side)
     144             : {
     145       86112 :   const Elem * current_remote_elem = _mesh.getMesh().elem_ptr(elem_id);
     146      172224 :   _current_remote_side = current_remote_elem->build_side_ptr(side);
     147       86112 :   _current_remote_side_volume = _assembly.elementVolume(_current_remote_side.get());
     148             : 
     149       86112 :   Order order = current_remote_elem->default_order();
     150       86112 :   unsigned int dim = _mesh.getMesh().mesh_dimension();
     151      172224 :   _current_remote_fe = FEBase::build(dim, FEType(order));
     152       86112 :   libMesh::QGauss qface(dim - 1, FEType(order).default_quadrature_order());
     153       86112 :   _current_remote_fe->attach_quadrature_rule(&qface);
     154             : 
     155       86112 :   _current_remote_JxW = &_current_remote_fe->get_JxW();
     156       86112 :   _current_remote_normals = &_current_remote_fe->get_normals();
     157       86112 :   _current_remote_xyz = &_current_remote_fe->get_xyz();
     158             : 
     159       86112 :   _current_remote_fe->reinit(current_remote_elem, side);
     160             : 
     161             :   // set _coord
     162       86112 :   unsigned int n_points = _current_remote_xyz->size();
     163       86112 :   unsigned int rz_radial_coord = _subproblem.getAxisymmetricRadialCoord();
     164       86112 :   _current_remote_coord.resize(n_points);
     165       86112 :   switch (_assembly.coordSystem())
     166             :   {
     167             :     case Moose::COORD_XYZ:
     168      362016 :       for (unsigned int qp = 0; qp < n_points; qp++)
     169      275904 :         _current_remote_coord[qp] = 1.;
     170             :       break;
     171             : 
     172             :     case Moose::COORD_RZ:
     173           0 :       for (unsigned int qp = 0; qp < n_points; qp++)
     174           0 :         _current_remote_coord[qp] = 2 * M_PI * (*_current_remote_xyz)[qp](rz_radial_coord);
     175             :       break;
     176             : 
     177             :     case Moose::COORD_RSPHERICAL:
     178           0 :       for (unsigned int qp = 0; qp < n_points; qp++)
     179           0 :         _current_remote_coord[qp] =
     180           0 :             4 * M_PI * (*_current_remote_xyz)[qp](0) * (*_current_remote_xyz)[qp](0);
     181             :       break;
     182             : 
     183           0 :     default:
     184           0 :       mooseError("Unknown coordinate system");
     185             :       break;
     186             :   }
     187       86112 : }