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