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 "PCNSFVKTDC.h"
      11             : #include "NS.h"
      12             : #include "SinglePhaseFluidProperties.h"
      13             : #include "Limiter.h"
      14             : #include "NavierStokesApp.h"
      15             : #include "Executioner.h"
      16             : 
      17             : using namespace Moose::FV;
      18             : 
      19             : registerMooseObject("NavierStokesApp", PCNSFVKTDC);
      20             : 
      21             : InputParameters
      22          76 : PCNSFVKTDC::validParams()
      23             : {
      24          76 :   InputParameters params = PCNSFVKT::validParams();
      25          76 :   params.addClassDescription("Computes the residual of advective term using finite volume method "
      26             :                              "using a deferred correction approach.");
      27         152 :   params.addParam<Real>(
      28         152 :       "ho_implicit_fraction", 0, "The fraction of the high order flux that should be implicit");
      29          76 :   return params;
      30           0 : }
      31             : 
      32          40 : PCNSFVKTDC::PCNSFVKTDC(const InputParameters & params)
      33             :   : PCNSFVKT(params),
      34          40 :     _upwind_limiter(Limiter<ADReal>::build(LimiterType::Upwind)),
      35          40 :     _old_upwind_fluxes(
      36          40 :         declareRestartableData<std::unordered_map<dof_id_type, Real>>("old_upwind_fluxes")),
      37          80 :     _old_ho_fluxes(declareRestartableData<std::unordered_map<dof_id_type, Real>>("old_ho_fluxes")),
      38          40 :     _current_upwind_fluxes(
      39          40 :         declareRestartableData<std::unordered_map<dof_id_type, Real>>("current_upwind_fluxes")),
      40          40 :     _current_ho_fluxes(
      41          40 :         declareRestartableData<std::unordered_map<dof_id_type, Real>>("current_ho_fluxes")),
      42         120 :     _ho_implicit_fraction(getParam<Real>("ho_implicit_fraction"))
      43             : {
      44          40 : }
      45             : 
      46             : void
      47         364 : PCNSFVKTDC::timestepSetup()
      48             : {
      49             :   // The SubProblem converged method is not restartable (it relies on data from libMesh) but the
      50             :   // Executioner converged method is, so we use that
      51         364 :   if (_app.getExecutioner()->lastSolveConverged())
      52             :   {
      53         360 :     _old_upwind_fluxes = _current_upwind_fluxes;
      54         360 :     _old_ho_fluxes = _current_ho_fluxes;
      55             :   }
      56         364 : }
      57             : 
      58             : void
      59        2548 : PCNSFVKTDC::residualSetup()
      60             : {
      61        2548 :   _current_upwind_fluxes.clear();
      62        2548 :   _current_ho_fluxes.clear();
      63        2548 : }
      64             : 
      65             : void
      66        1092 : PCNSFVKTDC::jacobianSetup()
      67             : {
      68        1092 :   residualSetup();
      69        1092 : }
      70             : 
      71             : Real
      72     1526952 : PCNSFVKTDC::getOldFlux(const bool upwind) const
      73             : {
      74     1526952 :   if (_t_step <= 1)
      75             :     return 0;
      76             : 
      77     1376760 :   const auto & flux_container = upwind ? _old_upwind_fluxes : _old_ho_fluxes;
      78     1376760 :   auto it = flux_container.find(_face_info->id());
      79             :   mooseAssert(it != flux_container.end(),
      80             :               "We should have saved an old flux for the current _face_info. Do you have mesh "
      81             :               "adaptivity on? Unfortunately we don't currently support that");
      82     1376760 :   return it->second;
      83             : }
      84             : 
      85             : ADReal
      86      763476 : PCNSFVKTDC::computeQpResidual()
      87             : {
      88             :   mooseAssert(!onBoundary(*_face_info), "We should never execute this object on a boundary");
      89             : 
      90      763476 :   const auto current_ho_flux = PCNSFVKT::computeQpResidual();
      91             : #ifndef NDEBUG
      92             :   auto pr =
      93             : #endif
      94      763476 :       _current_ho_fluxes.emplace(_face_info->id(), current_ho_flux.value());
      95             :   mooseAssert(pr.second,
      96             :               "Insertion should have happened. If it did not it means you are overwriting some "
      97             :               "other face's flux!");
      98             :   _limiter.swap(_upwind_limiter);
      99      763476 :   const auto current_upwind_flux = PCNSFVKT::computeQpResidual();
     100             : #ifndef NDEBUG
     101             :   pr =
     102             : #endif
     103      763476 :       _current_upwind_fluxes.emplace(_face_info->id(), current_upwind_flux.value());
     104             :   mooseAssert(pr.second,
     105             :               "Insertion should have happened. If it did not it means you are overwriting some "
     106             :               "other face's flux!");
     107             :   // Swap back
     108             :   _limiter.swap(_upwind_limiter);
     109             : 
     110      763476 :   const auto old_upwind_flux = getOldFlux(/*upwind=*/true);
     111      763476 :   const auto old_ho_flux = getOldFlux(/*upwind=*/false);
     112             : 
     113      763476 :   return _ho_implicit_fraction * current_ho_flux +
     114     1526952 :          (1. - _ho_implicit_fraction) * (old_ho_flux + current_upwind_flux - old_upwind_flux);
     115             : }