LCOV - code coverage report
Current view: top level - src/utils - CHTHandler.C (source / functions) Hit Total Coverage
Test: idaholab/moose navier_stokes: #32971 (54bef8) with base c6cf66 Lines: 204 224 91.1 %
Date: 2026-05-29 20:37:52 Functions: 11 11 100.0 %
Legend: Lines: hit not hit

          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             : // Moose includes
      11             : #include "CHTHandler.h"
      12             : #include "LinearFVFluxKernel.h"
      13             : #include "LinearFVDiffusion.h"
      14             : #include "LinearFVAnisotropicDiffusion.h"
      15             : #include "LinearFVBoundaryCondition.h"
      16             : #include "LinearFVCHTBCInterface.h"
      17             : #include "FEProblemBase.h"
      18             : 
      19             : namespace NS
      20             : {
      21             : namespace FV
      22             : {
      23             : 
      24             : InputParameters
      25        1344 : CHTHandler::validParams()
      26             : {
      27        1344 :   auto params = emptyInputParameters();
      28        2688 :   params.addParam<std::vector<BoundaryName>>(
      29             :       "cht_interfaces",
      30             :       {},
      31             :       "The interfaces where we would like to add conjugate heat transfer handling.");
      32             : 
      33        4032 :   params.addRangeCheckedParam<unsigned int>(
      34             :       "max_cht_fpi",
      35        2688 :       1,
      36             :       "max_cht_fpi >= 1",
      37             :       "Number of maximum fixed point iterations (FPI). Currently only applied to"
      38             :       " conjugate heat transfer simulations. The default value of 1 essentially keeps"
      39             :       " the FPI feature turned off. CHT iteration ends after this number of iteration even if the "
      40             :       "tolerance is not met.");
      41             : 
      42        4032 :   params.addRangeCheckedParam<Real>(
      43             :       "cht_heat_flux_tolerance",
      44        2688 :       1e-5,
      45             :       "cht_heat_flux_tolerance > 0 & cht_heat_flux_tolerance <= 1.0",
      46             :       "The relative tolerance for terminating conjugate heat transfer iteration before the maximum "
      47             :       "number of CHT iterations. Relative tolerance is ignore if the maximum number of CHT "
      48             :       "iterations is reached.");
      49             : 
      50        2688 :   params.addParam<std::vector<Real>>(
      51             :       "cht_fluid_temperature_relaxation",
      52             :       {},
      53             :       "The relaxation factors for the boundary temperature when being updated on the fluid side.");
      54        2688 :   params.addParam<std::vector<Real>>(
      55             :       "cht_solid_temperature_relaxation",
      56             :       {},
      57             :       "The relaxation factors for the boundary temperature when being updated on the solid side.");
      58        2688 :   params.addParam<std::vector<Real>>(
      59             :       "cht_fluid_flux_relaxation",
      60             :       {},
      61             :       "The relaxation factors for the boundary flux when being updated on the fluid side.");
      62        2688 :   params.addParam<std::vector<Real>>(
      63             :       "cht_solid_flux_relaxation",
      64             :       {},
      65             :       "The relaxation factors for the boundary flux when being updated on the solid side.");
      66             : 
      67        2688 :   params.addParamNamesToGroup(
      68             :       "cht_interfaces max_cht_fpi cht_heat_flux_tolerance cht_fluid_temperature_relaxation "
      69             :       "cht_solid_temperature_relaxation cht_fluid_flux_relaxation "
      70             :       "cht_solid_flux_relaxation",
      71             :       "Conjugate Heat Transfer");
      72             : 
      73        1344 :   return params;
      74           0 : }
      75             : 
      76         672 : CHTHandler::CHTHandler(const InputParameters & params)
      77             :   : MooseObject(params),
      78        1344 :     _problem(*getCheckedPointerParam<FEProblemBase *>(
      79             :         "_fe_problem_base", "This might happen if you don't have a mesh")),
      80         672 :     _mesh(_problem.mesh()),
      81        1344 :     _cht_boundary_names(getParam<std::vector<BoundaryName>>("cht_interfaces")),
      82         672 :     _cht_boundary_ids(_mesh.getBoundaryIDs(_cht_boundary_names)),
      83        1344 :     _max_cht_fpi(getParam<unsigned int>("max_cht_fpi")),
      84        2016 :     _cht_heat_flux_tolerance(getParam<Real>("cht_heat_flux_tolerance"))
      85             : {
      86        2016 :   if (isParamSetByUser("cht_interfaces") && !_cht_boundary_names.size())
      87           0 :     paramError("cht_interfaces", "You must declare at least one interface!");
      88         672 : }
      89             : 
      90             : void
      91          40 : CHTHandler::linkEnergySystems(SystemBase * solid_energy_system,
      92             :                               SystemBase * fluid_energy_system,
      93             :                               std::vector<SystemBase *> pm_radiation_systems)
      94             : {
      95          40 :   _energy_system = fluid_energy_system;
      96          40 :   _solid_energy_system = solid_energy_system;
      97          40 :   _pm_radiation_systems = pm_radiation_systems;
      98             : 
      99          40 :   if (!_energy_system || !_solid_energy_system)
     100           0 :     paramError("cht_interfaces",
     101             :                "You selected to do conjugate heat transfer treatment, but it needs two energy "
     102             :                "systems: a solid and a fluid. One of these systems is missing.");
     103          40 : }
     104             : 
     105             : void
     106          40 : CHTHandler::deduceCHTBoundaryCoupling()
     107             : {
     108          40 :   if (_solid_energy_system->nVariables() != 1)
     109           0 :     mooseError("We should have only one variable in the solid energy system: ",
     110           0 :                _energy_system->name(),
     111             :                "! Right now we have: ",
     112           0 :                Moose::stringify(_solid_energy_system->getVariableNames()));
     113          40 :   if (_energy_system->nVariables() != 1)
     114           0 :     mooseError("We should have only one variable in the fluid energy system: ",
     115           0 :                _energy_system->name(),
     116             :                "! Right now we have: ",
     117           0 :                Moose::stringify(_energy_system->getVariableNames()));
     118         120 :   const std::vector<std::string> solid_fluid({"solid", "fluid"});
     119             : 
     120             :   // We do some setup at the beginning to make sure the container sizes are good
     121             :   _cht_system_numbers =
     122          40 :       std::vector<unsigned int>({_solid_energy_system->number(), _energy_system->number()});
     123          40 :   _cht_conduction_kernels = std::vector<LinearFVFluxKernel *>({nullptr, nullptr});
     124          40 :   _cht_boundary_conditions.clear();
     125          40 :   _cht_boundary_conditions.resize(_cht_boundary_names.size(), {nullptr, nullptr});
     126             : 
     127             :   // Populate the PM radiation system numbers
     128          40 :   if (!_pm_radiation_systems.empty())
     129             :   {
     130          32 :     for (const auto sys_i : index_range(_pm_radiation_systems))
     131          16 :       _cht_pm_radiation_system_numbers.push_back(_pm_radiation_systems[sys_i]->number());
     132             : 
     133             :     // Reserve space for _cht_pm_radiation_kernels based on the size of
     134             :     // _cht_pm_radiation_system_numbers
     135          16 :     _cht_pm_radiation_kernels.reserve(_cht_pm_radiation_system_numbers.size());
     136             :     // Reserve space for pm radiation boundary conditions
     137          16 :     _cht_pm_radiation_boundary_conditions.clear();
     138          16 :     _cht_pm_radiation_boundary_conditions.resize(
     139             :         _cht_boundary_names.size(),
     140          32 :         std::vector<LinearFVBoundaryCondition *>(_cht_pm_radiation_system_numbers.size(), nullptr));
     141             :   }
     142             : 
     143             :   const auto flux_relaxation_param_names =
     144         120 :       std::vector<std::string>({"cht_solid_flux_relaxation", "cht_fluid_flux_relaxation"});
     145             :   const auto temperature_relaxation_param_names = std::vector<std::string>(
     146         120 :       {"cht_solid_temperature_relaxation", "cht_fluid_temperature_relaxation"});
     147          40 :   _cht_flux_relaxation_factor.clear();
     148          40 :   _cht_flux_relaxation_factor.resize(2, std::vector<Real>(_cht_boundary_names.size(), 1.0));
     149          40 :   _cht_temperature_relaxation_factor.clear();
     150          40 :   _cht_temperature_relaxation_factor.resize(2, std::vector<Real>(_cht_boundary_names.size(), 1.0));
     151             : 
     152         120 :   for (const auto region_index : index_range(solid_fluid))
     153             :   {
     154             :     // First thing, we fetch the relaxation parameter values
     155             :     const auto & flux_param_value =
     156             :         getParam<std::vector<Real>>(flux_relaxation_param_names[region_index]);
     157          80 :     if (flux_param_value.empty() || (flux_param_value.size() != _cht_boundary_names.size()))
     158           0 :       paramError(flux_relaxation_param_names[region_index],
     159             :                  "The number of relaxation factors is not the same as the number of interfaces!");
     160             : 
     161          80 :     _cht_flux_relaxation_factor[region_index] = flux_param_value;
     162             :     // We have to do the range check here because the intput parameter check errors if the vector is
     163             :     // empty
     164         160 :     for (const auto param : _cht_flux_relaxation_factor[region_index])
     165          80 :       if (param <= 0 || param > 1.0)
     166           0 :         paramError(flux_relaxation_param_names[region_index],
     167             :                    "The relaxation parameter should be between 0 and 1!");
     168             : 
     169             :     const auto & temperature_param_value =
     170             :         getParam<std::vector<Real>>(temperature_relaxation_param_names[region_index]);
     171          80 :     if (temperature_param_value.empty() ||
     172             :         (temperature_param_value.size() != _cht_boundary_names.size()))
     173           0 :       paramError(temperature_relaxation_param_names[region_index],
     174             :                  "The number of relaxation factors is not the same as the number of interfaces!");
     175             : 
     176          80 :     _cht_temperature_relaxation_factor[region_index] = temperature_param_value;
     177             :     // We have to do the range check here because the intput parameter check errors if the vector is
     178             :     // empty
     179         160 :     for (const auto param : _cht_temperature_relaxation_factor[region_index])
     180          80 :       if (param <= 0 || param > 1.0)
     181           0 :         paramError(temperature_relaxation_param_names[region_index],
     182             :                    "The relaxation parameter should be between 0 and 1!");
     183             : 
     184             :     // We then fetch the conduction kernels
     185             :     std::vector<LinearFVFluxKernel *> flux_kernels;
     186          80 :     _app.theWarehouse()
     187          80 :         .query()
     188          80 :         .template condition<AttribSystem>("LinearFVFluxKernel")
     189         160 :         .template condition<AttribVar>(0)
     190          80 :         .template condition<AttribSysNum>(_cht_system_numbers[region_index])
     191             :         .queryInto(flux_kernels);
     192             : 
     193             :     // We then fetch the radiation conduction kernels in the fluid region
     194          80 :     if (!_pm_radiation_systems.empty() && region_index == 1)
     195          32 :       for (const auto sys_i : index_range(_pm_radiation_systems))
     196             :       {
     197             :         // We then fetch the radiation conduction kernels
     198             :         std::vector<LinearFVFluxKernel *> radiation_kernels;
     199          32 :         _app.theWarehouse()
     200          16 :             .query()
     201          16 :             .template condition<AttribSystem>("LinearFVFluxKernel")
     202          32 :             .template condition<AttribVar>(0)
     203          16 :             .template condition<AttribSysNum>(_cht_pm_radiation_system_numbers[sys_i])
     204             :             .queryInto(radiation_kernels);
     205             : 
     206          16 :         if (radiation_kernels.size() > 1)
     207           0 :           mooseError(
     208             :               "We already have a kernel that describes the participating media radiation diffusion "
     209             :               "with the name: ",
     210           0 :               radiation_kernels[0]->name(),
     211             :               ". Make sure that you have only one conduction kernel.");
     212          16 :         else if (radiation_kernels.empty())
     213           0 :           mooseError("We did not find a diffusion kernel for the participating media radiation "
     214             :                      "diffusion to compute the "
     215             :                      "radiative heat flux. Please add a diffusion kernel.");
     216             :         else
     217          16 :           _cht_pm_radiation_kernels.push_back(radiation_kernels[0]);
     218          16 :       }
     219             : 
     220         184 :     for (auto kernel : flux_kernels)
     221             :     {
     222         104 :       auto check_diff = dynamic_cast<LinearFVDiffusion *>(kernel);
     223         104 :       auto check_aniso_diff = dynamic_cast<LinearFVAnisotropicDiffusion *>(kernel);
     224         104 :       if (_cht_conduction_kernels[region_index] && (check_diff || check_aniso_diff))
     225           0 :         mooseError("We already have a kernel that describes the heat conduction for the ",
     226             :                    solid_fluid[region_index],
     227             :                    " domain: ",
     228             :                    _cht_conduction_kernels[region_index]->name(),
     229             :                    " We found another one with the name: ",
     230             :                    (check_diff ? check_diff->name() : check_aniso_diff->name()),
     231             :                    " Make sure that you have only one conduction kernel on the ",
     232             :                    solid_fluid[region_index],
     233             :                    " side!");
     234             : 
     235         104 :       if (check_diff || check_aniso_diff)
     236          80 :         _cht_conduction_kernels[region_index] = kernel;
     237             :     }
     238             : 
     239             :     // Then we check the boundary conditions, to make sure at least there is something defined
     240             :     // from both sides
     241         160 :     for (const auto bd_index : index_range(_cht_boundary_names))
     242             :     {
     243             :       const auto & boundary_name = _cht_boundary_names[bd_index];
     244          80 :       const auto boundary_id = _cht_boundary_ids[bd_index];
     245             : 
     246             :       std::vector<LinearFVBoundaryCondition *> bcs;
     247          80 :       _problem.getMooseApp()
     248             :           .theWarehouse()
     249          80 :           .query()
     250          80 :           .template condition<AttribSystem>("LinearFVBoundaryCondition")
     251         160 :           .template condition<AttribVar>(0)
     252          80 :           .template condition<AttribSysNum>(_cht_system_numbers[region_index])
     253          80 :           .template condition<AttribBoundaries>(boundary_id)
     254             :           .queryInto(bcs);
     255             : 
     256             :       // We then fetch the radiation conduction bcs in the fluid region (i.e MarshakBC in P1)
     257          80 :       if (!_pm_radiation_systems.empty() && region_index == 1)
     258          32 :         for (const auto sys_i : index_range(_pm_radiation_systems))
     259             :         {
     260             :           std::vector<LinearFVBoundaryCondition *> rad_bcs;
     261          16 :           _app.theWarehouse()
     262          16 :               .query()
     263          16 :               .template condition<AttribSystem>("LinearFVBoundaryCondition")
     264          32 :               .template condition<AttribVar>(0)
     265          16 :               .template condition<AttribSysNum>(_cht_pm_radiation_system_numbers[sys_i])
     266          16 :               .template condition<AttribBoundaries>(boundary_id)
     267             :               .queryInto(rad_bcs);
     268             : 
     269          16 :           if (!rad_bcs.empty())
     270          16 :             _cht_pm_radiation_boundary_conditions[bd_index][sys_i] = rad_bcs[0];
     271             :           else
     272           0 :             mooseError("No LinearFVBoundaryCondition found for the given boundary or system.");
     273          16 :         }
     274             : 
     275          80 :       if (bcs.size() != 1)
     276           0 :         mooseError("We found multiple or no boundary conditions for solid energy on boundary ",
     277             :                    boundary_name,
     278             :                    " (ID: ",
     279             :                    boundary_id,
     280             :                    "). Make sure you define exactly one for conjugate heat transfer applications!");
     281          80 :       _cht_boundary_conditions[bd_index][region_index] = bcs[0];
     282             : 
     283          80 :       if (!dynamic_cast<LinearFVCHTBCInterface *>(_cht_boundary_conditions[bd_index][region_index]))
     284           0 :         mooseError("The selected boundary condition cannot be used with CHT problems! Make sure it "
     285             :                    "inherits from LinearFVCHTBCInterface!");
     286          80 :     }
     287          80 :   }
     288         240 : }
     289             : 
     290             : void
     291          40 : CHTHandler::setupConjugateHeatTransferContainers()
     292             : {
     293             :   // We already error in initialSetup if we have more variables
     294             :   const auto * fluid_variable =
     295          40 :       dynamic_cast<const MooseLinearVariableFVReal *>(&_energy_system->getVariable(0, 0));
     296             :   const auto * solid_variable =
     297          40 :       dynamic_cast<const MooseLinearVariableFVReal *>(&_solid_energy_system->getVariable(0, 0));
     298             : 
     299          40 :   _cht_face_info.clear();
     300          40 :   _cht_face_info.resize(_cht_boundary_ids.size());
     301          40 :   _boundary_heat_flux.clear();
     302          40 :   _boundary_temperature.clear();
     303          40 :   _integrated_boundary_heat_flux.clear();
     304             : 
     305          80 :   for (const auto bd_index : index_range(_cht_boundary_ids))
     306             :   {
     307          40 :     const auto bd_id = _cht_boundary_ids[bd_index];
     308             :     const auto & bd_name = _cht_boundary_names[bd_index];
     309             : 
     310             :     // We populate the face infos for every interface
     311             :     auto & bd_fi_container = _cht_face_info[bd_index];
     312       49734 :     for (auto & fi : _problem.mesh().faceInfo())
     313       49694 :       if (fi->boundaryIDs().count(bd_id))
     314         400 :         bd_fi_container.push_back(fi);
     315             : 
     316             :     // We do this because the coupling functors should be evaluated on both sides
     317             :     // of the interface and there are rigorous checks if the functors don't support a subdomain
     318             :     std::set<SubdomainID> combined_set;
     319          40 :     std::set_union(solid_variable->blockIDs().begin(),
     320          40 :                    solid_variable->blockIDs().end(),
     321          40 :                    fluid_variable->blockIDs().begin(),
     322          40 :                    fluid_variable->blockIDs().end(),
     323             :                    std::inserter(combined_set, combined_set.begin()));
     324             : 
     325             :     // We instantiate the coupling fuctors for heat flux and temperature
     326             :     FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> solid_bd_flux(
     327          40 :         _problem.mesh(), combined_set, "heat_flux_to_solid_" + bd_name);
     328             :     FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> fluid_bd_flux(
     329          40 :         _problem.mesh(), combined_set, "heat_flux_to_fluid_" + bd_name);
     330             : 
     331             :     _boundary_heat_flux.push_back(
     332         160 :         std::vector<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
     333          80 :             {std::move(solid_bd_flux), std::move(fluid_bd_flux)}));
     334             :     auto & flux_container = _boundary_heat_flux.back();
     335             : 
     336          80 :     _integrated_boundary_heat_flux.push_back(std::vector<Real>({0.0, 0.0}));
     337             : 
     338             :     FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> solid_bd_temperature(
     339          40 :         _problem.mesh(), combined_set, "interface_temperature_solid_" + bd_name);
     340             :     FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> fluid_bd_temperature(
     341          40 :         _problem.mesh(), combined_set, "interface_temperature_fluid_" + bd_name);
     342             : 
     343             :     _boundary_temperature.push_back(
     344         160 :         std::vector<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
     345          80 :             {std::move(solid_bd_temperature), std::move(fluid_bd_temperature)}));
     346             :     auto & temperature_container = _boundary_temperature.back();
     347             : 
     348             :     // Time to register the functors on all of the threads
     349          80 :     for (const auto tid : make_range(libMesh::n_threads()))
     350             :     {
     351          40 :       _problem.addFunctor("heat_flux_to_solid_" + bd_name, flux_container[NS::CHTSide::SOLID], tid);
     352          40 :       _problem.addFunctor("heat_flux_to_fluid_" + bd_name, flux_container[NS::CHTSide::FLUID], tid);
     353          40 :       _problem.addFunctor(
     354          40 :           "interface_temperature_solid_" + bd_name, temperature_container[NS::CHTSide::SOLID], tid);
     355          40 :       _problem.addFunctor(
     356          80 :           "interface_temperature_fluid_" + bd_name, temperature_container[NS::CHTSide::FLUID], tid);
     357             :     }
     358             : 
     359             :     // Initialize the containers, they will be filled with correct values soon.
     360             :     // Before any solve happens.
     361         120 :     for (const auto region_index : make_range(2))
     362         880 :       for (auto & fi : bd_fi_container)
     363             :       {
     364         800 :         flux_container[region_index][fi->id()] = 0.0;
     365         800 :         temperature_container[region_index][fi->id()] = 0.0;
     366             :       }
     367          40 :   }
     368         120 : }
     369             : 
     370             : void
     371          38 : CHTHandler::initializeCHTCouplingFields()
     372             : {
     373          76 :   for (const auto bd_index : index_range(_cht_boundary_ids))
     374             :   {
     375             :     const auto & bd_fi_container = _cht_face_info[bd_index];
     376             :     auto & temperature_container = _boundary_temperature[bd_index];
     377             : 
     378         114 :     for (const auto region_index : make_range(2))
     379             :     {
     380             :       // Can't be const considering we will update members from here
     381          76 :       auto bc = _cht_boundary_conditions[bd_index][region_index];
     382         872 :       for (const auto & fi : bd_fi_container)
     383             :       {
     384         796 :         bc->setupFaceData(fi, fi->faceType(std::make_pair(0, _cht_system_numbers[region_index])));
     385         796 :         temperature_container[1 - region_index][fi->id()] = bc->computeBoundaryValue();
     386             :       }
     387             :     }
     388             :   }
     389          38 : }
     390             : 
     391             : void
     392       14552 : CHTHandler::updateCHTBoundaryCouplingFields(const NS::CHTSide side)
     393             : {
     394             :   // Well we can just use the face that this enum casts into int very nicely
     395             :   // we can use it to get the index of the other side
     396       14552 :   const NS::CHTSide other_side = static_cast<NS::CHTSide>(1 - side);
     397             : 
     398       29104 :   for (const auto bd_index : index_range(_cht_boundary_ids))
     399             :   {
     400       14552 :     auto & other_bc = _cht_boundary_conditions[bd_index][other_side];
     401             :     auto & other_kernel = _cht_conduction_kernels[other_side];
     402             : 
     403             :     // We get the relaxation from the other side, so if we are fluid side we get the solid
     404             :     // relaxation
     405       14552 :     const auto temperature_relaxation = _cht_flux_relaxation_factor[other_side][bd_index];
     406       14552 :     const auto flux_relaxation = _cht_temperature_relaxation_factor[other_side][bd_index];
     407             : 
     408             :     // Fetching the right container here, if side is fluid we fetch "heat_flux_to_fluid"
     409       14552 :     auto & flux_container = _boundary_heat_flux[bd_index][side];
     410             :     // Fetching the other side's contaienr here, if side is fluid we fetch the solid temperature
     411             :     auto & temperature_container = _boundary_temperature[bd_index][other_side];
     412             :     // We will also update the integrated flux for output info
     413             :     auto & integrated_flux = _integrated_boundary_heat_flux[bd_index][side];
     414             :     // We are recomputing this so, time to zero this out
     415       14552 :     integrated_flux = 0.0;
     416             : 
     417             :     const auto & bd_fi_container = _cht_face_info[bd_index];
     418             : 
     419             :     // We enter the face loop to update the coupling fields
     420      181024 :     for (const auto & fi : bd_fi_container)
     421             :     {
     422      166472 :       other_kernel->setupFaceData(fi);
     423             :       // We will want the flux in W/m2 for the coupling so no face integral for now,
     424             :       // this can cause issues if we start using face area in the kernels
     425             :       // for more than just face integral multipliers.
     426             :       // Also, if we decide to not require overlapping meshes on the boundary
     427             :       // this will probably have to change.
     428      166472 :       other_kernel->setCurrentFaceArea(1.0);
     429      166472 :       other_bc->setupFaceData(fi, fi->faceType(std::make_pair(0, _cht_system_numbers[other_side])));
     430             : 
     431             :       // T_new = relaxation * T_boundary + (1-relaxation) * T_old
     432      332944 :       temperature_container[fi->id()] =
     433      166472 :           temperature_relaxation * other_bc->computeBoundaryValue() +
     434      166472 :           (1 - temperature_relaxation) * temperature_container[fi->id()];
     435             : 
     436             :       // Flux_new = relaxation * Flux_boundary + (1-relaxation) * Flux_old,
     437             :       // minus sign is due to the normal differences
     438             : 
     439             :       // Conductive flux
     440      166472 :       auto flux = other_kernel->computeBoundaryFlux(*other_bc);
     441             : 
     442             :       // If participating media radiation system exists we add the heat flux from the fluid
     443             :       // to the solid region.
     444      166472 :       if (!_pm_radiation_systems.empty() && side == NS::CHTSide::SOLID)
     445        4424 :         for (const auto sys_i : index_range(_pm_radiation_systems))
     446             :         {
     447        2212 :           _cht_pm_radiation_kernels[sys_i]->setupFaceData(fi);
     448        2212 :           _cht_pm_radiation_kernels[sys_i]->setCurrentFaceArea(1.0);
     449        2212 :           _cht_pm_radiation_boundary_conditions[bd_index][sys_i]->setupFaceData(
     450             :               fi, fi->faceType(std::make_pair(0, _cht_pm_radiation_system_numbers[sys_i])));
     451        2212 :           flux += _cht_pm_radiation_kernels[sys_i]->computeBoundaryFlux(
     452             :               *_cht_pm_radiation_boundary_conditions[bd_index][sys_i]);
     453             :         }
     454             : 
     455      332944 :       flux_container[fi->id()] =
     456      166472 :           flux_relaxation * flux + (1 - flux_relaxation) * flux_container[fi->id()];
     457             : 
     458             :       // We do the integral here
     459      166472 :       integrated_flux += flux * fi->faceArea() * fi->faceCoord();
     460             :     }
     461             :   }
     462       14552 : }
     463             : 
     464             : void
     465        7276 : CHTHandler::sumIntegratedFluxes()
     466             : {
     467       14552 :   for (const auto i : index_range(_integrated_boundary_heat_flux))
     468             :   {
     469             :     auto & integrated_fluxes = _integrated_boundary_heat_flux[i];
     470        7276 :     _problem.comm().sum(integrated_fluxes[NS::CHTSide::SOLID]);
     471        7276 :     _problem.comm().sum(integrated_fluxes[NS::CHTSide::FLUID]);
     472             :   }
     473        7276 : }
     474             : 
     475             : void
     476        7276 : CHTHandler::printIntegratedFluxes() const
     477             : {
     478       14552 :   for (const auto i : index_range(_integrated_boundary_heat_flux))
     479             :   {
     480             :     auto & integrated_fluxes = _integrated_boundary_heat_flux[i];
     481        7276 :     _console << " Iteration " << _fpi_it << " Boundary " << _cht_boundary_names[i]
     482        7276 :              << " flux on solid side " << integrated_fluxes[NS::CHTSide::SOLID]
     483        7276 :              << " flux on fluid side: " << integrated_fluxes[NS::CHTSide::FLUID] << std::endl;
     484             :   }
     485        7276 : }
     486             : 
     487             : void
     488        4572 : CHTHandler::resetIntegratedFluxes()
     489             : {
     490        9144 :   for (const auto i : index_range(_integrated_boundary_heat_flux))
     491        9144 :     _integrated_boundary_heat_flux[i] = std::vector<Real>({0.0, 0.0});
     492        4572 : }
     493             : 
     494             : bool
     495       56908 : CHTHandler::converged() const
     496             : {
     497       56908 :   if (_fpi_it >= _max_cht_fpi)
     498             :     return true;
     499             : 
     500       40582 :   for (const auto & boundary_flux : _integrated_boundary_heat_flux)
     501             :   {
     502       10378 :     const Real f1 = boundary_flux[0];
     503       10378 :     const Real f2 = boundary_flux[1];
     504             : 
     505             :     // Special case: both are zero at startup not converged yet
     506       10378 :     if (_fpi_it != 0 && (f1 == 0.0 && f2 == 0.0))
     507             :       return true;
     508             : 
     509             :     // These fluxes should be of opposite sign
     510       10378 :     const Real diff = std::abs(f1 + f2);
     511       10378 :     const Real denom = std::max({std::fabs(f1), std::fabs(f2), Real(1e-14)});
     512       10378 :     const Real rel_diff = diff / denom;
     513             : 
     514       10378 :     if (rel_diff >= _cht_heat_flux_tolerance)
     515             :       return false;
     516             :   }
     517             : 
     518       30204 :   return _fpi_it;
     519             : }
     520             : 
     521             : } // End FV namespace
     522             : } // End Moose namespace

Generated by: LCOV version 1.14