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 "ViewFactorBase.h"
      11             : #include "libmesh/quadrature.h"
      12             : 
      13             : #include <limits>
      14             : 
      15             : InputParameters
      16         559 : ViewFactorBase::validParams()
      17             : {
      18         559 :   InputParameters params = SideUserObject::validParams();
      19        1118 :   params.addParam<Real>("view_factor_tol",
      20        1118 :                         std::numeric_limits<Real>::max(),
      21             :                         "Tolerance for checking view factors. Default is to allow everything.");
      22        1118 :   params.addParam<bool>(
      23        1118 :       "print_view_factor_info", false, "Flag to print information about computed view factors.");
      24        1118 :   params.addParam<bool>("normalize_view_factor",
      25        1118 :                         true,
      26             :                         "Determines if view factors are normalized to sum to one (consistent with "
      27             :                         "their definition).");
      28         559 :   params.addClassDescription(
      29             :       "A base class for automatic computation of view factors between sidesets.");
      30         559 :   return params;
      31           0 : }
      32             : 
      33         295 : ViewFactorBase::ViewFactorBase(const InputParameters & parameters)
      34             :   : SideUserObject(parameters),
      35         590 :     _n_sides(boundaryIDs().size()),
      36         295 :     _areas(_n_sides),
      37         590 :     _view_factor_tol(getParam<Real>("view_factor_tol")),
      38         590 :     _normalize_view_factor(getParam<bool>("normalize_view_factor")),
      39         885 :     _print_view_factor_info(getParam<bool>("print_view_factor_info"))
      40             : {
      41             :   // sizing the view factor array
      42         295 :   _view_factors.resize(_n_sides);
      43        2293 :   for (auto & v : _view_factors)
      44        1998 :     v.resize(_n_sides);
      45             : 
      46             :   // set up the map from the side id to the local index & side name to local index
      47         590 :   std::vector<BoundaryName> boundary_names = getParam<std::vector<BoundaryName>>("boundary");
      48        2293 :   for (unsigned int j = 0; j < boundary_names.size(); ++j)
      49        1998 :     _side_name_index[boundary_names[j]] = j;
      50         295 : }
      51             : 
      52             : unsigned int
      53       11496 : ViewFactorBase::getSideNameIndex(std::string name) const
      54             : {
      55             :   auto it = _side_name_index.find(name);
      56       11496 :   if (it == _side_name_index.end())
      57           0 :     mooseError("Boundary ", name, " does not exist.");
      58       11496 :   return it->second;
      59             : }
      60             : 
      61             : Real
      62         284 : ViewFactorBase::getViewFactor(BoundaryID from_id, BoundaryID to_id) const
      63             : {
      64         284 :   auto from_name = _mesh.getBoundaryName(from_id);
      65         284 :   auto to_name = _mesh.getBoundaryName(to_id);
      66             : 
      67         852 :   return getViewFactor(from_name, to_name);
      68             : }
      69             : 
      70             : Real
      71      537145 : ViewFactorBase::getViewFactor(BoundaryName from_name, BoundaryName to_name) const
      72             : {
      73      537145 :   auto from = _side_name_index.find(from_name);
      74      537145 :   auto to = _side_name_index.find(to_name);
      75      537145 :   if (from == _side_name_index.end())
      76           0 :     mooseError("Boundary id ",
      77           0 :                _mesh.getBoundaryID(from_name),
      78             :                " with name ",
      79             :                from_name,
      80             :                " not listed in boundary parameter.");
      81             : 
      82      537145 :   if (to == _side_name_index.end())
      83           0 :     mooseError("Boundary id ",
      84           0 :                _mesh.getBoundaryID(to_name),
      85             :                " with name ",
      86             :                to_name,
      87             :                " not listed in boundary parameter.");
      88             : 
      89      537145 :   return _view_factors[from->second][to->second];
      90             : }
      91             : 
      92             : void
      93         194 : ViewFactorBase::finalize()
      94             : {
      95             :   // do some communication before finalizing view_factors
      96        1510 :   for (unsigned int i = 0; i < _n_sides; ++i)
      97        1316 :     gatherSum(_view_factors[i]);
      98             : 
      99         194 :   finalizeViewFactor();
     100         194 :   checkAndNormalizeViewFactor();
     101         194 : }
     102             : 
     103             : void
     104          31 : ViewFactorBase::threadJoin(const UserObject & y)
     105             : {
     106             :   const auto & pps = static_cast<const ViewFactorBase &>(y);
     107         236 :   for (unsigned int i = 0; i < _n_sides; ++i)
     108             :   {
     109        1638 :     for (unsigned int j = 0; j < _n_sides; ++j)
     110        1433 :       _view_factors[i][j] += pps._view_factors[i][j];
     111             :   }
     112          31 :   threadJoinViewFactor(y);
     113          31 : }
     114             : 
     115             : Real
     116       19368 : ViewFactorBase::devReciprocity(unsigned int i, unsigned int j) const
     117             : {
     118       19368 :   return _view_factors[i][j] - _areas[j] / _areas[i] * _view_factors[j][i];
     119             : }
     120             : 
     121             : Real
     122         388 : ViewFactorBase::maxDevReciprocity() const
     123             : {
     124             :   Real v = 0;
     125        3020 :   for (unsigned int i = 0; i < _n_sides; ++i)
     126       10776 :     for (unsigned int j = i + 1; j < _n_sides; ++j)
     127             :     {
     128        8144 :       Real s = std::abs(devReciprocity(i, j));
     129        8144 :       if (s > v)
     130             :         v = s;
     131             :     }
     132         388 :   return v;
     133             : }
     134             : 
     135             : Real
     136        3450 : ViewFactorBase::viewFactorRowSum(unsigned int i) const
     137             : {
     138             :   Real s = 0;
     139       28800 :   for (unsigned int to = 0; to < _n_sides; ++to)
     140       25350 :     s += _view_factors[i][to];
     141        3450 :   return s;
     142             : }
     143             : 
     144             : Real
     145         388 : ViewFactorBase::maxDevRowSum() const
     146             : {
     147             :   Real v = 0;
     148        3020 :   for (unsigned int i = 0; i < _n_sides; ++i)
     149             :   {
     150        2632 :     Real s = std::abs(1 - viewFactorRowSum(i));
     151        2632 :     if (s > v)
     152             :       v = s;
     153             :   }
     154         388 :   return v;
     155             : }
     156             : 
     157             : void
     158         194 : ViewFactorBase::checkAndNormalizeViewFactor()
     159             : {
     160             :   // collect statistics
     161         194 :   Real max_sum_deviation = maxDevRowSum();
     162         194 :   Real max_reciprocity_deviation = maxDevReciprocity();
     163         194 :   Real max_correction = std::numeric_limits<Real>::lowest();
     164         194 :   Real min_correction = std::numeric_limits<Real>::max();
     165             : 
     166         194 :   if (_print_view_factor_info)
     167           0 :     _console << "\nSum of all view factors from side i to all other faces before normalization.\n"
     168           0 :              << std::endl;
     169             : 
     170             :   // check view factors
     171         194 :   if (_print_view_factor_info)
     172           0 :     for (unsigned int from = 0; from < _n_sides; ++from)
     173           0 :       _console << "View factors from sideset " << boundaryNames()[from] << " sum to "
     174           0 :                << viewFactorRowSum(from) << std::endl;
     175             : 
     176         194 :   if (max_sum_deviation > _view_factor_tol)
     177           0 :     mooseError("Maximum deviation of view factor row sum is ",
     178             :                max_sum_deviation,
     179             :                " exceeding the set tolerance of ",
     180           0 :                _view_factor_tol);
     181             : 
     182             :   // normalize view factors
     183         194 :   if (_normalize_view_factor)
     184             :   {
     185         110 :     if (_print_view_factor_info)
     186           0 :       _console << "\nNormalizing view factors.\n" << std::endl;
     187             : 
     188             :     // allocate space
     189         110 :     DenseVector<Real> rhs(_n_sides);
     190         110 :     DenseVector<Real> lmults(_n_sides);
     191         110 :     DenseMatrix<Real> matrix(_n_sides, _n_sides);
     192             : 
     193             :     // equations for the Lagrange multiplier
     194         928 :     for (unsigned int i = 0; i < _n_sides; ++i)
     195             :     {
     196             :       // set the right hand side
     197         818 :       rhs(i) = 1 - viewFactorRowSum(i);
     198             : 
     199             :       // contribution from the delta_ii element
     200         818 :       matrix(i, i) = -0.5;
     201             : 
     202             :       // contributions from the upper diagonal
     203        3624 :       for (unsigned int j = i + 1; j < _n_sides; ++j)
     204             :       {
     205        2806 :         Real ar = _areas[i] / _areas[j];
     206        2806 :         Real f = 2 * (1 + ar * ar);
     207        2806 :         matrix(i, i) += -1 / f;
     208        2806 :         matrix(i, j) += -1 / f * ar;
     209        2806 :         rhs(i) += ar * ar / (1 + ar * ar) * devReciprocity(i, j);
     210             :       }
     211             : 
     212             :       // contributions from the lower diagonal
     213        3624 :       for (unsigned int j = 0; j < i; ++j)
     214             :       {
     215        2806 :         Real ar = _areas[j] / _areas[i];
     216        2806 :         Real f = 2 * (1 + ar * ar);
     217        2806 :         matrix(i, i) += -1 / f * ar * ar;
     218        2806 :         matrix(i, j) += -1 / f * ar;
     219        2806 :         rhs(i) -= ar * devReciprocity(j, i) * (1 - ar * ar / (1 + ar * ar));
     220             :       }
     221             :     }
     222             : 
     223             :     // solve the linear system
     224         110 :     matrix.lu_solve(rhs, lmults);
     225             : 
     226             :     // perform corrections but store view factors to temporary array
     227             :     // because we cannot modify _view_factors as it's used in this calc
     228         110 :     std::vector<std::vector<Real>> vf_temp = _view_factors;
     229         928 :     for (unsigned int i = 0; i < _n_sides; ++i)
     230        7248 :       for (unsigned int j = 0; j < _n_sides; ++j)
     231             :       {
     232             :         Real correction;
     233        6430 :         if (i == j)
     234         818 :           correction = -0.5 * lmults(i);
     235             :         else
     236             :         {
     237        5612 :           Real ar = _areas[i] / _areas[j];
     238        5612 :           Real f = 2 * (1 + ar * ar);
     239        5612 :           correction = -(lmults(i) + lmults(j) * ar + 2 * ar * ar * devReciprocity(i, j)) / f;
     240             :         }
     241             : 
     242             :         // update the temporary view factor
     243        6430 :         vf_temp[i][j] += correction;
     244             : 
     245        6430 :         if (correction < min_correction)
     246         608 :           min_correction = correction;
     247        6430 :         if (correction > max_correction)
     248         587 :           max_correction = correction;
     249             :       }
     250         110 :     _view_factors = vf_temp;
     251         110 :   }
     252             : 
     253         194 :   if (_print_view_factor_info)
     254             :   {
     255           0 :     _console << "\nFinal view factors.\n" << std::endl;
     256           0 :     for (unsigned int from = 0; from < _n_sides; ++from)
     257             :     {
     258             :       std::string from_name;
     259           0 :       for (auto pair : _side_name_index)
     260           0 :         if (pair.second == from)
     261             :           from_name = pair.first;
     262           0 :       auto from_id = _mesh.getBoundaryID(from_name);
     263             : 
     264           0 :       for (unsigned int to = 0; to < _n_sides; ++to)
     265             :       {
     266             :         std::string to_name;
     267           0 :         for (auto pair : _side_name_index)
     268           0 :           if (pair.second == to)
     269             :             to_name = pair.first;
     270           0 :         auto to_id = _mesh.getBoundaryID(to_name);
     271           0 :         _console << from_name << " (" << from_id << ") -> " << to_name << " (" << to_id
     272           0 :                  << ") = " << _view_factors[from][to] << std::endl;
     273             :       }
     274             :     }
     275             :   }
     276             : 
     277             :   // check sum of view factors and maximum deviation on reciprocity
     278         194 :   Real max_sum_deviation_after_normalization = maxDevRowSum();
     279         194 :   Real max_reciprocity_deviation_after_normalization = maxDevReciprocity();
     280             : 
     281             :   // print summary
     282         194 :   _console << std::endl;
     283         225 :   _console << COLOR_CYAN << "Summary of the view factor computation"
     284         194 :            << "\n"
     285         225 :            << COLOR_DEFAULT << std::endl;
     286         194 :   if (_normalize_view_factor)
     287         110 :     _console << "Normalization is performed." << std::endl;
     288             :   else
     289          84 :     _console << "Normalization is skipped." << std::endl;
     290         194 :   _console << std::setw(60) << std::left << "Number of patches: " << _n_sides << std::endl;
     291         194 :   _console << std::setw(60) << std::left
     292         194 :            << "Max deviation sum = 1 before normalization: " << max_sum_deviation << std::endl;
     293         194 :   _console << std::setw(60) << std::left
     294         194 :            << "Max deviation from reciprocity before normalization: " << max_reciprocity_deviation
     295         194 :            << std::endl;
     296         194 :   if (_normalize_view_factor)
     297             :   {
     298         110 :     _console << std::setw(60) << std::left << "Maximum correction: " << max_correction << std::endl;
     299         110 :     _console << std::setw(60) << std::left << "Minimum correction: " << min_correction << std::endl;
     300         110 :     _console << std::setw(60) << "Max deviation sum = 1 after normalization: "
     301         110 :              << max_sum_deviation_after_normalization << std::endl;
     302         110 :     _console << std::setw(60) << std::left << "Max deviation from reciprocity after normalization: "
     303         110 :              << max_reciprocity_deviation_after_normalization << std::endl;
     304         110 :     _console << std::endl;
     305             :   }
     306         194 : }
     307             : 
     308             : unsigned int
     309           0 : ViewFactorBase::indexHelper(unsigned int i, unsigned int j) const
     310             : {
     311             :   mooseAssert(i <= j, "indexHelper requires i <= j");
     312           0 :   if (i == j)
     313           0 :     return _n_sides + i;
     314           0 :   unsigned int pos = 2 * _n_sides;
     315           0 :   for (unsigned int l = 0; l < _n_sides; ++l)
     316           0 :     for (unsigned int m = l + 1; m < _n_sides; ++m)
     317             :     {
     318           0 :       if (l == i && m == j)
     319           0 :         return pos;
     320             :       else
     321           0 :         ++pos;
     322             :     }
     323           0 :   mooseError("Should never get here");
     324             :   return 0;
     325             : }