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 "DiscreteNucleationTimeStep.h"
      11             : #include "DiscreteNucleationInserterBase.h"
      12             : #include "MooseUtils.h"
      13             : 
      14             : registerMooseObject("PhaseFieldApp", DiscreteNucleationTimeStep);
      15             : 
      16             : InputParameters
      17          38 : DiscreteNucleationTimeStep::validParams()
      18             : {
      19          38 :   InputParameters params = GeneralPostprocessor::validParams();
      20          38 :   params.addClassDescription(
      21             :       "Return a time step limit for nucleation event to be used by IterationAdaptiveDT");
      22          76 :   params.addRequiredParam<Real>("dt_max",
      23             :                                 "Time step to cut back to at the start of a nucleation event");
      24         114 :   params.addRangeCheckedParam<Real>(
      25             :       "p2nucleus",
      26          76 :       0.01,
      27             :       "p2nucleus > 0 & p2nucleus < 1",
      28             :       "Maximum probability for more than one nucleus to appear during a time "
      29             :       "step. This will limit the time step based on the total nucleation rate for "
      30             :       "the domain to make sure the probability for two or more nuclei to appear "
      31             :       "is always below the chosen number.");
      32          76 :   params.addRequiredParam<UserObjectName>("inserter", "DiscreteNucleationInserter user object");
      33          38 :   return params;
      34           0 : }
      35             : 
      36          19 : DiscreteNucleationTimeStep::DiscreteNucleationTimeStep(const InputParameters & parameters)
      37             :   : GeneralPostprocessor(parameters),
      38          19 :     _inserter(getUserObject<DiscreteNucleationInserterBase>("inserter")),
      39          38 :     _dt_nucleation(getParam<Real>("dt_max")),
      40          19 :     _changes_made(_inserter.getInsertionsAndDeletions()),
      41          38 :     _rate(_inserter.getRate())
      42             : {
      43             :   //
      44             :   // We do a bisection search because math is hard
      45             :   // (i.e. probability function is not analytically invertible)
      46             :   //
      47             : 
      48             :   // this is the target value
      49          38 :   const Real p2n = getParam<Real>("p2nucleus");
      50             : 
      51             :   // initial guess
      52          19 :   _max_lambda = 0.1;
      53             :   Real upper_bound = _max_lambda;
      54             :   Real lower_bound = 0.0;
      55             : 
      56             :   // At this point we do not know a proper upper bound for a search interval
      57             :   // so we grow our initial guess until we find it (we only allow for a finite
      58             :   // number of iterations)
      59          52 :   for (unsigned int i = 0; i < 100; ++i)
      60             :   {
      61          52 :     const Real p_upper = 1.0 - (1.0 + upper_bound) * std::exp(-upper_bound);
      62             : 
      63             :     // we have found a lambda value that results in a p > p2n, this is the upper end of the interval
      64          52 :     if (p_upper > p2n)
      65             :       break;
      66             : 
      67             :     // upper_bound is actually below our target lambda, set it as the new lower
      68             :     // bound and double the upper_bound value
      69             :     lower_bound = upper_bound;
      70          33 :     upper_bound *= 2.0;
      71             :   }
      72             : 
      73             :   // now that we have an upper and a lower interval bounds we can do a proper bisection
      74         636 :   for (unsigned int i = 0; i < 100; ++i)
      75             :   {
      76             :     // pick the middle of the current interval
      77         636 :     _max_lambda = (upper_bound - lower_bound) / 2.0 + lower_bound;
      78             : 
      79             :     // calculate new probability for 2 or more nuclei to appear
      80         636 :     const Real p = 1.0 - (1.0 + _max_lambda) * std::exp(-_max_lambda);
      81             : 
      82             :     // quit if we zeroed in on the target
      83         636 :     if (MooseUtils::absoluteFuzzyEqual(p, p2n))
      84             :       break;
      85             : 
      86             :     // otherwise adjust interval bounds
      87         617 :     else if (p > p2n)
      88             :       upper_bound = _max_lambda;
      89             :     else
      90             :       lower_bound = _max_lambda;
      91             :   }
      92          19 : }
      93             : 
      94             : PostprocessorValue
      95         161 : DiscreteNucleationTimeStep::getValue() const
      96             : {
      97             :   // check if a nucleus insertion has occurred...
      98         161 :   if (_changes_made.first > 0)
      99             :     // and cut back time step for nucleus formation
     100          40 :     return _dt_nucleation;
     101             : 
     102             :   // otherwise check the total nucleation rate in the domain...
     103         121 :   if (_rate == 0.0)
     104             :     // ...and return no limit on the time step if the rate is zero...
     105             :     return std::numeric_limits<Real>::max();
     106             :   else
     107             :     // ...or return the maximum time step that satisfies the bound on the 2+ nuclei probability
     108         121 :     return _max_lambda / _rate;
     109             : }