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