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 "DiscreteNucleationInserter.h"
      11             : #include "libmesh/parallel_algebra.h"
      12             : #include "SystemBase.h"
      13             : 
      14             : #include "libmesh/quadrature.h"
      15             : 
      16             : registerMooseObject("PhaseFieldApp", DiscreteNucleationInserter);
      17             : 
      18             : InputParameters
      19         276 : DiscreteNucleationInserter::validParams()
      20             : {
      21         276 :   InputParameters params = DiscreteNucleationInserterBase::validParams();
      22         276 :   params.addClassDescription("Manages the list of currently active nucleation sites and adds new "
      23             :                              "sites according to a given probability function.");
      24         552 :   params.addRequiredParam<MaterialPropertyName>(
      25             :       "probability", "Probability density for inserting a discrete nucleus");
      26         552 :   params.addRequiredParam<Real>("hold_time", "Time to keep each nucleus active");
      27         552 :   params.addParam<MaterialPropertyName>("radius",
      28             :                                         "r_crit",
      29             :                                         "variable radius material property name, supply a value if "
      30             :                                         "radius is constant in the simulation");
      31         552 :   params.addParam<bool>("time_dependent_statistics",
      32         552 :                         true,
      33             :                         "flag if time-dependent or time-independent statistics are used");
      34             : 
      35         276 :   return params;
      36           0 : }
      37             : 
      38         146 : DiscreteNucleationInserter::DiscreteNucleationInserter(const InputParameters & parameters)
      39             :   : DiscreteNucleationInserterBase(parameters),
      40         146 :     _probability(getMaterialProperty<Real>("probability")),
      41         292 :     _hold_time(getParam<Real>("hold_time")),
      42         292 :     _local_nucleus_list(declareRestartableData<NucleusList>("local_nucleus_list", 0)),
      43         292 :     _local_radius(getMaterialProperty<Real>("radius")),
      44         438 :     _time_dep_stats(getParam<bool>("time_dependent_statistics"))
      45             : {
      46         146 : }
      47             : 
      48             : void
      49        1018 : DiscreteNucleationInserter::initialize()
      50             : {
      51             :   // clear insertion and deletion counter
      52             :   _changes_made = {0, 0};
      53             : 
      54             :   // expire entries from the local nucleus list (if the current time step converged)
      55        1018 :   if (_fe_problem.converged(_sys.number()))
      56             :   {
      57             :     unsigned int i = 0;
      58        1815 :     while (i < _local_nucleus_list.size())
      59             :     {
      60         797 :       if (_local_nucleus_list[i].time <= _fe_problem.time())
      61             :       {
      62             :         // remove entry (by replacing with last element and shrinking size by one)
      63         151 :         _local_nucleus_list[i] = _local_nucleus_list.back();
      64             :         _local_nucleus_list.pop_back();
      65         151 :         _changes_made.second++;
      66             :       }
      67             :       else
      68         646 :         ++i;
      69             :     }
      70             :   }
      71             : 
      72             :   // we reassemble this list at every time step
      73        1018 :   _global_nucleus_list.clear();
      74             : 
      75             :   // clear total nucleation rate
      76        1018 :   _nucleation_rate = 0.0;
      77        1018 : }
      78             : 
      79             : void
      80      969368 : DiscreteNucleationInserter::execute()
      81             : {
      82             :   // check each qp for potential nucleation
      83             :   // TODO: we might as well place the nuclei at random positions within the element...
      84     4774840 :   for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
      85             :   {
      86     3805472 :     const Real rate = _probability[qp] * _JxW[qp] * _coord[qp];
      87     3805472 :     _nucleation_rate += rate;
      88             : 
      89     3805472 :     const Real random = getRandomReal();
      90             : 
      91             :     // branch the operation for using time-dependent statistics or
      92             :     // time-independent probability (e.g., recrystallization fraction)
      93             :     // If time-dependent, `rate` refers to a probability rate density
      94             :     // If time-independent, `rate` refers to a probability density
      95     3805472 :     if (!_time_dep_stats)
      96             :     {
      97             :       // if using time-independent statistics, this would be more like a nucleation fraction
      98           0 :       if (random < rate)
      99           0 :         addNucleus(qp);
     100             :     }
     101             :     else
     102             :     {
     103             :       // We check the random number against the inverse of the zero probability.
     104             :       // for performance reasons we do a quick check against the linearized form of
     105             :       // that probability, which is always strictly larger than the actual probability.
     106             :       // The expression below should short circuit and the expensive exponential
     107             :       // should rarely get evaluated
     108     3805472 :       if (random < rate * _fe_problem.dt() && random < (1.0 - std::exp(-rate * _fe_problem.dt())))
     109         315 :         addNucleus(qp);
     110             :     }
     111             :   }
     112      969368 : }
     113             : 
     114             : void
     115         130 : DiscreteNucleationInserter::threadJoin(const UserObject & y)
     116             : {
     117             :   // combine _local_nucleus_list entries from all threads on the current process
     118             :   const auto & uo = static_cast<const DiscreteNucleationInserter &>(y);
     119         260 :   _global_nucleus_list.insert(
     120         130 :       _global_nucleus_list.end(), uo._local_nucleus_list.begin(), uo._local_nucleus_list.end());
     121             : 
     122             :   // sum up insertion and deletion counts
     123         130 :   _changes_made.first += uo._changes_made.first;
     124         130 :   _changes_made.second += uo._changes_made.second;
     125             : 
     126             :   // integrate total nucleation rate
     127         130 :   _nucleation_rate += uo._nucleation_rate;
     128         130 : }
     129             : 
     130             : void
     131         888 : DiscreteNucleationInserter::finalize()
     132             : {
     133             :   // add the _local_nucleus_list of thread zero
     134        1776 :   _global_nucleus_list.insert(
     135         888 :       _global_nucleus_list.end(), _local_nucleus_list.begin(), _local_nucleus_list.end());
     136             : 
     137             :   /**
     138             :    * Pack the _global_nucleus_list into a simple vector of Real.
     139             :    * libMesh's allgather does not portably work on the original
     140             :    * _global_nucleus_list data structure!
     141             :    */
     142         888 :   std::vector<Real> comm_buffer(_global_nucleus_list.size() * 5);
     143        1849 :   for (unsigned i = 0; i < _global_nucleus_list.size(); ++i)
     144             :   {
     145         961 :     comm_buffer[i * 5 + 0] = _global_nucleus_list[i].time;
     146         961 :     comm_buffer[i * 5 + 1] = _global_nucleus_list[i].center(0);
     147         961 :     comm_buffer[i * 5 + 2] = _global_nucleus_list[i].center(1);
     148         961 :     comm_buffer[i * 5 + 3] = _global_nucleus_list[i].center(2);
     149         961 :     comm_buffer[i * 5 + 4] = _global_nucleus_list[i].radius;
     150             :   }
     151             : 
     152             :   // combine _global_nucleus_lists from all MPI ranks
     153         888 :   _communicator.allgather(comm_buffer);
     154             : 
     155             :   // unpack the gathered _global_nucleus_list
     156         888 :   unsigned int n = comm_buffer.size() / 5;
     157             :   mooseAssert(comm_buffer.size() % 5 == 0,
     158             :               "Communication buffer has an unexpected size (not divisible by 5)");
     159         888 :   _global_nucleus_list.resize(n);
     160             : 
     161        2864 :   for (unsigned i = 0; i < n; ++i)
     162             :   {
     163        1976 :     _global_nucleus_list[i].time = comm_buffer[i * 5 + 0];
     164        1976 :     _global_nucleus_list[i].center(0) = comm_buffer[i * 5 + 1];
     165        1976 :     _global_nucleus_list[i].center(1) = comm_buffer[i * 5 + 2];
     166        1976 :     _global_nucleus_list[i].center(2) = comm_buffer[i * 5 + 3];
     167        1976 :     _global_nucleus_list[i].radius = comm_buffer[i * 5 + 4];
     168             :   }
     169             : 
     170             :   // get the global number of changes (i.e. changes to _global_nucleus_list)
     171         888 :   gatherSum(_changes_made.first);
     172         888 :   gatherSum(_changes_made.second);
     173             : 
     174             :   // gather the total nucleation rate
     175         888 :   gatherSum(_nucleation_rate);
     176             : 
     177         888 :   _update_required = _changes_made.first > 0 || _changes_made.second > 0;
     178         888 : }
     179             : 
     180             : void
     181         315 : DiscreteNucleationInserter::addNucleus(unsigned int & qp)
     182             : {
     183             :   NucleusLocation new_nucleus;
     184         315 :   new_nucleus.time = _fe_problem.time() + _hold_time;
     185         315 :   new_nucleus.center = _q_point[qp];
     186         315 :   new_nucleus.radius = _local_radius[qp];
     187             : 
     188         315 :   _local_nucleus_list.push_back(new_nucleus);
     189         315 :   _changes_made.first++;
     190         315 : }