Alamo
PhaseFieldMicrostructure.cpp
Go to the documentation of this file.
1 
2 #include <eigen3/Eigen/Eigenvalues>
3 
4 #include <cmath>
5 
6 #include <AMReX_SPACE.H>
7 
8 #include "PhaseFieldMicrostructure.H"
9 #include "Integrator/Base/Mechanics.H"
10 #include "BC/Constant.H"
11 #include "Set/Set.H"
12 #include "Util/Util.H"
13 #include "IC/Random.H"
14 #include "IC/Trig.H"
15 #include "IC/Sphere.H"
16 #include "IC/Expression.H"
17 #include "Model/Interface/GB/SH.H"
18 #include "Numeric/Stencil.H"
19 #include "Solver/Nonlocal/Linear.H"
20 #include "Solver/Nonlocal/Newton.H"
21 #include "IC/Trig.H"
22 
23 #include "Model/Solid/Affine/Cubic.H"
24 #include "Model/Solid/Affine/Hexagonal.H"
25 
26 namespace Integrator
27 {
28 template<class model_type>
29 void PhaseFieldMicrostructure<model_type>::Advance(int lev, Set::Scalar time, Set::Scalar dt)
30 {
31  BL_PROFILE("PhaseFieldMicrostructure::Advance");
32  Base::Mechanics<model_type>::Advance(lev, time, dt);
33  /// TODO Make this optional
34  //if (lev != max_level) return;
35  //std::swap(eta_old_mf[lev], eta_new_mf[lev]);
36  const Set::Scalar* DX = this->geom[lev].CellSize();
37 
38 
39  Model::Interface::GB::SH gbmodel(0.0, 0.0, anisotropy.sigma0, anisotropy.sigma1);
40 
41  for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
42  {
43  amrex::Box bx = mfi.tilebox();
44  //if (m_type == MechanicsBase<model_type>::Type::Static)
45  //bx.grow(number_of_ghost_cells-1);
46  //bx = bx & domain;
47  Set::Patch<Set::Scalar> eta = eta_mf.Patch(lev,mfi);
48  Set::Patch<Set::Scalar> driving_force = driving_force_mf.Patch(lev,mfi);
49  Set::Patch<Set::Scalar> driving_force_threshold = driving_force_threshold_mf.Patch(lev,mfi);// = (*driving_force_mf[lev]).array(mfi);
50 
51  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
52  {
53  for (int m = 0; m < number_of_grains; m++)
54  {
55  driving_force(i, j, k, m) = 0.0;
56  if (pf.threshold.on) driving_force_threshold(i, j, k, m) = 0.0;
57  Set::Scalar kappa = NAN, mu = NAN;
58 
59  //
60  // BOUNDARY TERM and SECOND ORDER REGULARIZATION
61  //
62 
63  Set::Vector Deta = Numeric::Gradient(eta, i, j, k, m, DX);
64  Set::Scalar normgrad = Deta.lpNorm<2>();
65  if (normgrad < 1E-4)
66  continue; // This ought to speed things up.
67 
68  Set::Matrix DDeta = Numeric::Hessian(eta, i, j, k, m, DX);
69  Set::Scalar laplacian = DDeta.trace();
70 
71  if (!anisotropy.on || time < anisotropy.tstart)
72  {
73  kappa = pf.l_gb * 0.75 * pf.sigma0;
74  mu = 0.75 * (1.0 / 0.23) * pf.sigma0 / pf.l_gb;
75  if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += -kappa * laplacian;
76  else driving_force(i, j, k, m) += -kappa * laplacian;
77  }
78  else
79  {
80  Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Full> DDDDEta = Numeric::DoubleHessian<AMREX_SPACEDIM>(eta, i, j, k, m, DX);
81  auto anisotropic_df = boundary->DrivingForce(Deta, DDeta, DDDDEta);
82  if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += pf.l_gb * 0.75 * std::get<0>(anisotropic_df);
83  else driving_force(i, j, k, m) += pf.l_gb * 0.75 * std::get<0>(anisotropic_df);
84  if (std::isnan(std::get<0>(anisotropic_df))) Util::Abort(INFO);
85  if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += anisotropy.beta * std::get<1>(anisotropic_df);
86  else driving_force(i, j, k, m) += anisotropy.beta * std::get<1>(anisotropic_df);
87  if (std::isnan(std::get<1>(anisotropic_df))) Util::Abort(INFO);
88  mu = 0.75 * (1.0 / 0.23) * boundary->W(Deta) / pf.l_gb;
89  }
90 
91  //
92  // CHEMICAL POTENTIAL
93  //
94 
95  Set::Scalar sum_of_squares = 0.;
96  for (int n = 0; n < number_of_grains; n++)
97  {
98  if (m == n)
99  continue;
100  sum_of_squares += eta(i, j, k, n) * eta(i, j, k, n);
101  }
102  if (pf.threshold.chempot)
103  driving_force_threshold(i, j, k, m) += mu * (eta(i, j, k, m) * eta(i, j, k, m) - 1.0 + 2.0 * pf.gamma * sum_of_squares) * eta(i, j, k, m);
104  else
105  driving_force(i, j, k, m) += mu * (eta(i, j, k, m) * eta(i, j, k, m) - 1.0 + 2.0 * pf.gamma * sum_of_squares) * eta(i, j, k, m);
106 
107  //
108  // LAGRANGE MULTIPLIER
109  //
110  if (lagrange.on && m == 0 && time > lagrange.tstart)
111  {
112  if (pf.threshold.lagrange)
113  driving_force_threshold(i, j, k, m) += lagrange.lambda * (volume - lagrange.vol0);
114  else
115  driving_force(i, j, k, m) += lagrange.lambda * (volume - lagrange.vol0);
116  }
117 
118  //
119  // SYNTHETIC DRIVING FORCE
120  //
121  if (sdf.on && time > sdf.tstart)
122  {
123  if (pf.threshold.sdf)
124  driving_force_threshold(i, j, k, m) += sdf.val[m](time);
125  else
126  driving_force(i, j, k, m) += sdf.val[m](time);
127  }
128  }
129  });
130 
131  //
132  // ELASTIC DRIVING FORCE
133  //
134  if (pf.elastic_df)
135  {
136  amrex::Array4<const Set::Matrix> const& sigma = (*this->stress_mf[lev]).array(mfi);
137 
138  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
139  {
140  for (int m = 0; m < number_of_grains; m++)
141  {
142  Set::Matrix F0avg = Set::Matrix::Zero();
143 
144  for (int n = 0; n < number_of_grains; n++)
145  {
146  F0avg += eta(i, j, k, n) * mechanics.model[n].F0;
147  }
148 
149  Set::Matrix sig = Numeric::Interpolate::NodeToCellAverage(sigma, i, j, k, 0);
150 
151  //Set::Matrix dF0deta = mechanics.model[m].F0;//(etasum * elastic.model[m].F0 - F0avg) / (etasum * etasum);
152  Set::Matrix dF0deta = Set::Matrix::Zero();
153 
154  for (int n = 0; n < number_of_grains; n++)
155  {
156  if (n == m) continue;
157  Set::Scalar normsq = eta(i, j, k, m) * eta(i, j, k, m) + eta(i, j, k, n) * eta(i, j, k, n);
158  dF0deta += (2.0 * eta(i, j, k, m) * eta(i, j, k, n) * eta(i, j, k, n) * (mechanics.model[m].F0 - mechanics.model[n].F0))
159  / normsq / normsq;
160  }
161 
162  Set::Scalar tmpdf = (dF0deta.transpose() * sig).trace();
163 
164  if (pf.threshold.mechanics)
165  driving_force_threshold(i, j, k, m) -= pf.elastic_mult * tmpdf;
166  else
167  driving_force(i, j, k, m) -= pf.elastic_mult * tmpdf;
168  }
169  });
170  }
171 
172 
173  //
174  // Update eta
175  // (if NOT using anisotropic kinetics)
176  //
177 
178  if (!anisotropic_kinetics.on || time < anisotropic_kinetics.tstart)
179  {
180  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
181  {
182  for (int m = 0; m < number_of_grains; m++)
183  {
184  if (pf.threshold.on)
185  {
186  if (driving_force_threshold(i, j, k, m) > pf.threshold.value)
187  {
188  if (pf.threshold.type == ThresholdType::Continuous)
189  eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m) - pf.threshold.value);
190  else if (pf.threshold.type == ThresholdType::Chop)
191  eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m));
192  }
193  else if (driving_force_threshold(i, j, k, m) < -pf.threshold.value)
194  {
195  if (pf.threshold.type == ThresholdType::Continuous)
196  eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m) + pf.threshold.value);
197  else if (pf.threshold.type == ThresholdType::Chop)
198  eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m));
199  }
200  }
201 
202  eta(i, j, k, m) -= pf.L * dt * driving_force(i, j, k, m);
203  }
204  });
205  }
206  }
207 
208  //
209  // Update eta
210  // (if we ARE using anisotropic kinetics)
211  //
212 
213  if (anisotropic_kinetics.on && time >= anisotropic_kinetics.tstart)
214  {
215  for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
216  {
217  amrex::Box bx = mfi.tilebox();
218  amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
219  amrex::Array4<Set::Scalar> const& L = (*anisotropic_kinetics.L_mf[lev]).array(mfi);
220  amrex::Array4<Set::Scalar> const& threshold = (*anisotropic_kinetics.threshold_mf[lev]).array(mfi);
221 
222  for (int m = 0; m < number_of_grains; m++)
223  {
224  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
225  {
226  Set::Vector Deta = Numeric::Gradient(eta, i, j, k, m, DX);
227  Set::Scalar theta = atan2(Deta(1), Deta(0));
228  L(i, j, k, m) = (4. / 3.) * anisotropic_kinetics.mobility(theta) / pf.l_gb;
229  threshold(i, j, k, m) = anisotropic_kinetics.threshold(theta);
230  });
231  }
232  }
233 
234  for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
235  {
236  amrex::Box bx = mfi.tilebox();
237  Set::Patch<Set::Scalar> eta = eta_mf.Patch(lev,mfi);
238  Set::Patch<const Set::Scalar> driving_force = driving_force_mf.Patch(lev,mfi);
239  Set::Patch<const Set::Scalar> driving_force_threshold = driving_force_threshold_mf.Patch(lev,mfi);
240  Set::Patch<const Set::Scalar> L = anisotropic_kinetics.L_mf.Patch(lev,mfi);
241  Set::Patch<const Set::Scalar> threshold = anisotropic_kinetics.threshold_mf.Patch(lev,mfi);
242 
243  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
244  {
245  for (int m = 0; m < number_of_grains; m++)
246  {
247  if (pf.threshold.on)
248  {
249  if (driving_force_threshold(i, j, k, m) > threshold(i,j,k,m))
250  {
251  if (pf.threshold.type == ThresholdType::Continuous)
252  eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m) - threshold(i,j,k,m));
253  else if (pf.threshold.type == ThresholdType::Chop)
254  eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m));
255  }
256  else if (driving_force_threshold(i, j, k, m) < -pf.threshold.value)
257  {
258  if (pf.threshold.type == ThresholdType::Continuous)
259  eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m) + threshold(i,j,k,m));
260  else if (pf.threshold.type == ThresholdType::Chop)
261  eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m));
262  }
263  }
264 
265  eta(i, j, k, m) -= L(i,j,k,m) * dt * driving_force(i, j, k, m);
266  }
267  });
268  }
269  }
270 }
271 
272 template<class model_type>
273 void PhaseFieldMicrostructure<model_type>::Initialize(int lev)
274 {
275  BL_PROFILE("PhaseFieldMicrostructure::Initialize");
276  Base::Mechanics<model_type>::Initialize(lev);
277  ic->Initialize(lev, eta_mf);
278 }
279 
280 template<class model_type>
281 void PhaseFieldMicrostructure<model_type>::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
282 {
283  BL_PROFILE("PhaseFieldMicrostructure::TagCellsForRefinement");
284  Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
285  const amrex::Real* DX = this->geom[lev].CellSize();
286  const Set::Vector dx(DX);
287  const Set::Scalar dxnorm = dx.lpNorm<2>();
288 
289  for (amrex::MFIter mfi(*eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
290  {
291  const amrex::Box& bx = mfi.tilebox();
292  amrex::Array4<const amrex::Real> const& etanew = (*eta_mf[lev]).array(mfi);
293  amrex::Array4<char> const& tags = a_tags.array(mfi);
294 
295  for (int n = 0; n < number_of_grains; n++)
296  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
297  Set::Vector grad = Numeric::Gradient(etanew, i, j, k, n, DX);
298 
299  if (dxnorm * grad.lpNorm<2>() > ref_threshold)
300  tags(i, j, k) = amrex::TagBox::SET;
301  });
302  }
303 }
304 
305 template<class model_type>
306 void PhaseFieldMicrostructure<model_type>::TimeStepComplete(Set::Scalar /*time*/, int /*iter*/)
307 {
308  // TODO: remove this function, it is no longer needed.
309 }
310 
311 template<class model_type>
312 void PhaseFieldMicrostructure<model_type>::UpdateModel(int a_step, Set::Scalar /*a_time*/)
313 {
314  BL_PROFILE("PhaseFieldMicrostructure::UpdateModel");
315  if (a_step % this->m_interval) return;
316 
317  for (int lev = 0; lev <= this->finest_level; ++lev)
318  {
319  amrex::Box domain = this->geom[lev].Domain();
320  domain.convert(amrex::IntVect::TheNodeVector());
321 
322  eta_mf[lev]->FillBoundary();
323 
324  for (MFIter mfi(*this->model_mf[lev], false); mfi.isValid(); ++mfi)
325  {
326  amrex::Box bx = mfi.grownnodaltilebox() & domain;
327 
328  amrex::Array4<model_type> const& model = this->model_mf[lev]->array(mfi);
329  amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
330 
331  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
332  {
333  std::vector<Set::Scalar> etas(number_of_grains);
334  for (int n = 0; n < number_of_grains; n++)
335  etas[n] = Numeric::Interpolate::CellToNodeAverage(eta, i, j, k, n);
336  model(i, j, k) = model_type::Combine(mechanics.model, etas, 1);
337  });
338  }
339 
340  Util::RealFillBoundary(*this->model_mf[lev], this->geom[lev]);
341  }
342 
343 }
344 
345 template<class model_type>
346 void PhaseFieldMicrostructure<model_type>::TimeStepBegin(Set::Scalar time, int iter)
347 {
348  BL_PROFILE("PhaseFieldMicrostructure::TimeStepBegin");
349  Base::Mechanics<model_type>::TimeStepBegin(time, iter);
350 
351  if (anisotropy.on && time >= anisotropy.tstart)
352  {
353  this->SetTimestep(anisotropy.timestep);
354  if (anisotropy.elastic_int > 0)
355  if (iter % anisotropy.elastic_int) return;
356  }
357 }
358 
359 template<class model_type>
360 void PhaseFieldMicrostructure<model_type>::Integrate(int amrlev, Set::Scalar time, int step,
361  const amrex::MFIter& mfi, const amrex::Box& box)
362 {
363  BL_PROFILE("PhaseFieldMicrostructure::Integrate");
364  Base::Mechanics<model_type>::Integrate(amrlev, time, step, mfi, box);
365 
366  Model::Interface::GB::SH gbmodel(0.0, 0.0, anisotropy.sigma0, anisotropy.sigma1);
367  const amrex::Real* DX = this->geom[amrlev].CellSize();
368  Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
369 
370  amrex::Array4<amrex::Real> const& eta = (*eta_mf[amrlev]).array(mfi);
371  amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
372 #if AMREX_SPACEDIM == 2
373  auto sten = Numeric::GetStencil(i, j, k, box);
374 #endif
375 
376  volume += eta(i, j, k, 0) * dv;
377 
378  Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
379  Set::Scalar normgrad = grad.lpNorm<2>();
380 
381  if (normgrad > 1E-8)
382  {
383  Set::Vector normal = grad / normgrad;
384 
385  Set::Scalar da = normgrad * dv;
386  area += da;
387 
388  if (!anisotropy.on || time < anisotropy.tstart)
389  {
390  gbenergy += pf.sigma0 * da;
391  Set::Scalar k = 0.75 * pf.sigma0 * pf.l_gb;
392  realgbenergy += 0.5 * k * normgrad * normgrad * dv;
393  regenergy = 0.0;
394  }
395  else
396  {
397 #if AMREX_SPACEDIM == 2
398  Set::Scalar theta = atan2(grad(1), grad(0));
399  Set::Scalar sigma = boundary->W(theta);
400  gbenergy += sigma * da;
401 
402  Set::Scalar k = 0.75 * sigma * pf.l_gb;
403  realgbenergy += 0.5 * k * normgrad * normgrad * dv;
404 
405  Set::Matrix DDeta = Numeric::Hessian(eta, i, j, k, 0, DX, sten);
406  Set::Vector tangent(normal[1], -normal[0]);
407  Set::Scalar k2 = (DDeta * tangent).dot(tangent);
408  regenergy += 0.5 * anisotropy.beta * k2 * k2;
409 #elif AMREX_SPACEDIM == 3
410  gbenergy += gbmodel.W(normal) * da;
411 #endif
412  }
413  }
414  });
415 }
416 
417 template class PhaseFieldMicrostructure<Model::Solid::Affine::Cubic>;
418 template class PhaseFieldMicrostructure<Model::Solid::Affine::Hexagonal>;
419 
420 
421 } // namespace Integrator
Util::Initialize
void Initialize()
Definition: Util.cpp:126
SH.H
Numeric::Interpolate::CellToNodeAverage
static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:1277
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T >> &a_mf, const amrex::Geometry &)
Definition: Util.H:319
Integrator::Base::Mechanics
Definition: Mechanics.H:23
Util.H
Sphere.H
Expression.H
Model::Interface::GB::SH::W
Set::Scalar W(const Set::Scalar a_theta, const Set::Scalar a_phi) const
Definition: SH.H:60
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
Constant.H
Trig.H
Numeric::Gradient
AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:656
Integrator::PhaseFieldMicrostructure
Microstructure evolution with grain boundary anisotropy.
Definition: PhaseFieldMicrostructure.H:57
Cubic.H
Numeric::Interpolate::NodeToCellAverage
static AMREX_FORCE_INLINE T NodeToCellAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m)
Definition: Stencil.H:1315
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition: Base.H:23
Mechanics.H
PhaseFieldMicrostructure.H
Model::Interface::GB::SH
Definition: SH.H:35
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:166
Integrator
Collection of numerical integrator objects.
Definition: AllenCahn.H:40
Random.H
Newton.H
Set::Matrix4
Definition: Base.H:198
Set.H
Numeric::Hessian
AMREX_FORCE_INLINE Set::Matrix Hessian(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition: Stencil.H:930
Integrator::PhaseFieldMicrostructure::Advance
void Advance(int lev, Real time, Real dt) override
Evolve phase field in time.
Definition: PhaseFieldMicrostructure.cpp:29
Integrator::Base::Mechanics::Advance
void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
Definition: Mechanics.H:232
Set::Patch
amrex::Array4< T > const & Patch
Definition: Set.H:88
Hexagonal.H
INFO
#define INFO
Definition: Util.H:20
Numeric::GetStencil
static AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition: Stencil.H:20
Linear.H
Generated on Sat Feb 8 2025 00:07:39 for Alamo by   doxygen 1.8.17