www.mooseframework.org
DiscreteNucleationInserter.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
11 #include "libmesh/parallel_algebra.h"
12 #include "SystemBase.h"
13 
14 #include "libmesh/quadrature.h"
15 
17 
20 {
22  params.addClassDescription("Manages the list of currently active nucleation sites and adds new "
23  "sites according to a given probability function.");
24  params.addRequiredParam<MaterialPropertyName>(
25  "probability", "Probability density for inserting a discrete nucleus");
26  params.addRequiredParam<Real>("hold_time", "Time to keep each nucleus active");
27  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  params.addParam<bool>("time_dependent_statistics",
32  true,
33  "flag if time-dependent or time-independent statistics are used");
34 
35  return params;
36 }
37 
39  : DiscreteNucleationInserterBase(parameters),
40  _probability(getMaterialProperty<Real>("probability")),
41  _hold_time(getParam<Real>("hold_time")),
42  _local_nucleus_list(declareRestartableData<NucleusList>("local_nucleus_list", 0)),
43  _local_radius(getMaterialProperty<Real>("radius")),
44  _time_dep_stats(getParam<bool>("time_dependent_statistics"))
45 {
46 }
47 
48 void
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)
56  {
57  unsigned int i = 0;
58  while (i < _local_nucleus_list.size())
59  {
60  if (_local_nucleus_list[i].time <= _fe_problem.time())
61  {
62  // remove entry (by replacing with last element and shrinking size by one)
64  _local_nucleus_list.pop_back();
65  _changes_made.second++;
66  }
67  else
68  ++i;
69  }
70  }
71 
72  // we reassemble this list at every time step
73  _global_nucleus_list.clear();
74 
75  // clear total nucleation rate
76  _nucleation_rate = 0.0;
77 }
78 
79 void
81 {
82  // check each qp for potential nucleation
83  // TODO: we might as well place the nuclei at random positions within the element...
84  for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
85  {
86  const Real rate = _probability[qp] * _JxW[qp] * _coord[qp];
87  _nucleation_rate += rate;
88 
89  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  if (!_time_dep_stats)
96  {
97  // if using time-independent statistics, this would be more like a nucleation fraction
98  if (random < rate)
99  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  if (random < rate * _fe_problem.dt() && random < (1.0 - std::exp(-rate * _fe_problem.dt())))
109  addNucleus(qp);
110  }
111  }
112 }
113 
114 void
116 {
117  // combine _local_nucleus_list entries from all threads on the current process
118  const auto & uo = static_cast<const DiscreteNucleationInserter &>(y);
119  _global_nucleus_list.insert(
120  _global_nucleus_list.end(), uo._local_nucleus_list.begin(), uo._local_nucleus_list.end());
121 
122  // sum up insertion and deletion counts
123  _changes_made.first += uo._changes_made.first;
124  _changes_made.second += uo._changes_made.second;
125 
126  // integrate total nucleation rate
127  _nucleation_rate += uo._nucleation_rate;
128 }
129 
130 void
132 {
133  // add the _local_nucleus_list of thread zero
134  _global_nucleus_list.insert(
136 
142  std::vector<Real> comm_buffer(_global_nucleus_list.size() * 5);
143  for (unsigned i = 0; i < _global_nucleus_list.size(); ++i)
144  {
145  comm_buffer[i * 5 + 0] = _global_nucleus_list[i].time;
146  comm_buffer[i * 5 + 1] = _global_nucleus_list[i].center(0);
147  comm_buffer[i * 5 + 2] = _global_nucleus_list[i].center(1);
148  comm_buffer[i * 5 + 3] = _global_nucleus_list[i].center(2);
149  comm_buffer[i * 5 + 4] = _global_nucleus_list[i].radius;
150  }
151 
152  // combine _global_nucleus_lists from all MPI ranks
153  _communicator.allgather(comm_buffer);
154 
155  // unpack the gathered _global_nucleus_list
156  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  _global_nucleus_list.resize(n);
160 
161  for (unsigned i = 0; i < n; ++i)
162  {
163  _global_nucleus_list[i].time = comm_buffer[i * 5 + 0];
164  _global_nucleus_list[i].center(0) = comm_buffer[i * 5 + 1];
165  _global_nucleus_list[i].center(1) = comm_buffer[i * 5 + 2];
166  _global_nucleus_list[i].center(2) = comm_buffer[i * 5 + 3];
167  _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  gatherSum(_changes_made.first);
172  gatherSum(_changes_made.second);
173 
174  // gather the total nucleation rate
176 
177  _update_required = _changes_made.first > 0 || _changes_made.second > 0;
178 }
179 
180 void
182 {
183  NucleusLocation new_nucleus;
184  new_nucleus.time = _fe_problem.time() + _hold_time;
185  new_nucleus.center = _q_point[qp];
186  new_nucleus.radius = _local_radius[qp];
187 
188  _local_nucleus_list.push_back(new_nucleus);
189  _changes_made.first++;
190 }
NucleusChanges _changes_made
count the number of nucleus insertions and deletions
void allgather(const T &send_data, std::vector< T, A > &recv_data) const
registerMooseObject("PhaseFieldApp", DiscreteNucleationInserter)
virtual Real & time() const
const MooseArray< Point > & _q_point
virtual bool converged(const unsigned int nl_sys_num)
const MooseArray< Real > & _coord
virtual void addNucleus(unsigned int &qp)
Adds a nucleus to the list containing nucleus information.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
const MaterialProperty< Real > & _local_radius
store the local nucleus radius
const MaterialProperty< Real > & _probability
Nucleation rate density (should be a material property implementing nucleation theory) ...
std::vector< NucleusLocation > NucleusList
Every MPI task should keep a full list of nuclei (in case they cross domains with their finite radii)...
A nucleus has an expiration time, a location, and a size.
const std::vector< double > y
const Parallel::Communicator & _communicator
static InputParameters validParams()
void addRequiredParam(const std::string &name, const std::string &doc_string)
This UserObject manages the insertion and expiration of nuclei in the simulation domain it manages a ...
This UserObject manages the insertion and expiration of nuclei in the simulation domain it manages a ...
DiscreteNucleationInserter(const InputParameters &parameters)
Real _nucleation_rate
Total nucleation rate.
void gatherSum(T &value)
virtual void threadJoin(const UserObject &y)
NucleusList & _local_nucleus_list
the local nucleus list of nuclei centered in the domain of the current processor
SystemBase & _sys
unsigned int number() const
bool _update_required
is a map update required
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real _hold_time
Duration of time each nucleus is kept active after insertion.
const QBase *const & _qrule
FEProblemBase & _fe_problem
const MooseArray< Real > & _JxW
NucleusList & _global_nucleus_list
the global list of all nuclei over all processors
void addClassDescription(const std::string &doc_string)
virtual Real & dt() const
Real getRandomReal() const
const bool _time_dep_stats
indicates whether time-dependent statistics are used or not