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 "INSExplicitTimestepSelector.h"
      11             : 
      12             : #include "libmesh/quadrature.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", INSExplicitTimestepSelector);
      15             : 
      16             : InputParameters
      17           0 : INSExplicitTimestepSelector::validParams()
      18             : {
      19           0 :   InputParameters params = ElementPostprocessor::validParams();
      20             : 
      21           0 :   params.addClassDescription("Postprocessor that computes the minimum value of h_min/|u|, where "
      22             :                              "|u| is coupled in as an aux variable.");
      23             : 
      24             :   // Coupled variables
      25           0 :   params.addRequiredCoupledVar("vel_mag", "Velocity magnitude");
      26             : 
      27             :   // Required parameters
      28           0 :   params.addRequiredParam<Real>("beta",
      29             :                                 "0 < beta < 1, choose some fraction of the limiting timestep size");
      30             : 
      31             :   // Optional parameters
      32           0 :   params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
      33           0 :   params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
      34             : 
      35           0 :   return params;
      36           0 : }
      37             : 
      38           0 : INSExplicitTimestepSelector::INSExplicitTimestepSelector(const InputParameters & parameters)
      39             :   : ElementPostprocessor(parameters),
      40           0 :     _vel_mag(coupledValue("vel_mag")),
      41             : 
      42             :     // Other parameters
      43           0 :     _beta(getParam<Real>("beta")),
      44             : 
      45             :     // Material properties
      46           0 :     _mu(getMaterialProperty<Real>("mu_name")),
      47           0 :     _rho(getMaterialProperty<Real>("rho_name"))
      48             : {
      49           0 : }
      50             : 
      51           0 : INSExplicitTimestepSelector::~INSExplicitTimestepSelector() {}
      52             : 
      53             : void
      54           0 : INSExplicitTimestepSelector::initialize()
      55             : {
      56           0 :   _value = std::numeric_limits<Real>::max();
      57           0 : }
      58             : 
      59             : void
      60           0 : INSExplicitTimestepSelector::execute()
      61             : {
      62           0 :   Real h_min = _current_elem->hmin();
      63             : 
      64             :   // The space dimension plays a role in the diffusive dt limit.  The more
      65             :   // space dimensions there are, the smaller this limit is.
      66           0 :   Real dim = static_cast<Real>(_current_elem->dim());
      67             : 
      68           0 :   for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
      69             :   {
      70             :     // Don't divide by zero...
      71           0 :     Real vel_mag = std::max(_vel_mag[qp], std::numeric_limits<Real>::epsilon());
      72             : 
      73             :     // For explicit Euler, we always have to satisfy the Courant condition for stability.
      74           0 :     Real courant_limit_dt = h_min / vel_mag;
      75             : 
      76             :     // But we also have to obey the diffusive time restriction,
      77             :     // dt <= 1/(2*nu)*(1/h1^2 + 1/h2^2 + 1/h3^2)^(-1) <=
      78             :     //    <= h_min^2 / n_dim / (2*nu)
      79           0 :     Real diffusive_limit_dt = 0.5 * h_min * h_min / (_mu[qp] / _rho[qp]) / dim;
      80             : 
      81             :     // And the "combined" condition, dt <= 2*nu/|u|^2
      82           0 :     Real combined_limit_dt = 2. * (_mu[qp] / _rho[qp]) / vel_mag / vel_mag;
      83             : 
      84             :     // // Debugging:
      85             :     // Moose::out << "courant_limit_dt   = " << courant_limit_dt   << "\n"
      86             :     //           << "diffusive_limit_dt = " << diffusive_limit_dt << "\n"
      87             :     //           << "combined_limit_dt  = " << combined_limit_dt
      88             :     //           << std::endl;
      89             : 
      90           0 :     _value = std::min(
      91           0 :         _value,
      92           0 :         _beta * std::min(std::min(courant_limit_dt, diffusive_limit_dt), combined_limit_dt));
      93             :   }
      94           0 : }
      95             : 
      96             : Real
      97           0 : INSExplicitTimestepSelector::getValue() const
      98             : {
      99           0 :   return _value;
     100             : }
     101             : 
     102             : void
     103           0 : INSExplicitTimestepSelector::finalize()
     104             : {
     105           0 :   gatherMin(_value);
     106           0 : }
     107             : 
     108             : void
     109           0 : INSExplicitTimestepSelector::threadJoin(const UserObject & uo)
     110             : {
     111             :   const auto & pps = static_cast<const INSExplicitTimestepSelector &>(uo);
     112           0 :   _value = std::min(_value, pps._value);
     113           0 : }