Line data Source code
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 324 : void PhaseFieldMicrostructure<model_type>::Advance(int lev, Set::Scalar time, Set::Scalar dt)
31 : {
32 : BL_PROFILE("PhaseFieldMicrostructure::Advance");
33 324 : 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 324 : const Set::Scalar* DX = this->geom[lev].CellSize();
38 :
39 :
40 648 : Model::Interface::GB::SH gbmodel(0.0, 0.0, anisotropy.sigma0, anisotropy.sigma1);
41 :
42 30986 : for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
43 : {
44 30662 : 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 30662 : Set::Patch<Set::Scalar> eta = eta_mf.Patch(lev,mfi);
49 30662 : Set::Patch<Set::Scalar> driving_force = driving_force_mf.Patch(lev,mfi);
50 30662 : Set::Patch<Set::Scalar> driving_force_threshold = driving_force_threshold_mf.Patch(lev,mfi);// = (*driving_force_mf[lev]).array(mfi);
51 :
52 2029062 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
53 : {
54 10005700 : for (int m = 0; m < number_of_grains; m++)
55 : {
56 8007300 : driving_force(i, j, k, m) = 0.0;
57 8007300 : if (pf.threshold.on) driving_force_threshold(i, j, k, m) = 0.0;
58 8007300 : Set::Scalar kappa = NAN, mu = NAN;
59 :
60 : //
61 : // BOUNDARY TERM and SECOND ORDER REGULARIZATION
62 : //
63 :
64 8007300 : Set::Vector Deta = Numeric::Gradient(eta, i, j, k, m, DX);
65 8007300 : Set::Scalar normgrad = Deta.lpNorm<2>();
66 8007300 : if (normgrad < 1E-4)
67 5766470 : continue; // This ought to speed things up.
68 :
69 2240830 : Set::Matrix DDeta = Numeric::Hessian(eta, i, j, k, m, DX);
70 2240830 : Set::Scalar laplacian = DDeta.trace();
71 :
72 2240830 : if (!anisotropy.on || time < anisotropy.tstart)
73 : {
74 877969 : kappa = pf.l_gb * 0.75 * pf.sigma0;
75 877969 : mu = 0.75 * (1.0 / 0.23) * pf.sigma0 / pf.l_gb;
76 877969 : if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += -kappa * laplacian;
77 1755940 : else driving_force(i, j, k, m) += -kappa * laplacian;
78 : }
79 : else
80 : {
81 1362860 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Full> DDDDEta = Numeric::DoubleHessian<AMREX_SPACEDIM>(eta, i, j, k, m, DX);
82 1362860 : auto anisotropic_df = boundary->DrivingForce(Deta, DDeta, DDDDEta);
83 1362860 : if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += pf.l_gb * 0.75 * std::get<0>(anisotropic_df);
84 2725720 : else driving_force(i, j, k, m) += pf.l_gb * 0.75 * std::get<0>(anisotropic_df);
85 1362860 : if (std::isnan(std::get<0>(anisotropic_df))) Util::Abort(INFO);
86 1362860 : if (pf.threshold.boundary) driving_force_threshold(i, j, k, m) += anisotropy.beta * std::get<1>(anisotropic_df);
87 2725720 : else driving_force(i, j, k, m) += anisotropy.beta * std::get<1>(anisotropic_df);
88 1362860 : if (std::isnan(std::get<1>(anisotropic_df))) Util::Abort(INFO);
89 1362860 : mu = 0.75 * (1.0 / 0.23) * boundary->W(Deta) / pf.l_gb;
90 : }
91 :
92 : //
93 : // CHEMICAL POTENTIAL
94 : //
95 :
96 2240830 : Set::Scalar sum_of_squares = 0.;
97 13398500 : for (int n = 0; n < number_of_grains; n++)
98 : {
99 11157600 : if (m == n)
100 2240830 : continue;
101 26750400 : sum_of_squares += eta(i, j, k, n) * eta(i, j, k, n);
102 : }
103 2240830 : if (pf.threshold.chempot)
104 0 : 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 11204100 : 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 2240830 : if (lagrange.on && m == 0 && time > lagrange.tstart)
112 : {
113 0 : if (pf.threshold.lagrange)
114 0 : driving_force_threshold(i, j, k, m) += lagrange.lambda * (volume - lagrange.vol0);
115 : else
116 0 : driving_force(i, j, k, m) += lagrange.lambda * (volume - lagrange.vol0);
117 : }
118 :
119 : //
120 : // SYNTHETIC DRIVING FORCE
121 : //
122 2240830 : if (sdf.on && time > sdf.tstart)
123 : {
124 0 : if (pf.threshold.sdf)
125 0 : driving_force_threshold(i, j, k, m) += sdf.val[m](time);
126 : else
127 0 : driving_force(i, j, k, m) += sdf.val[m](time);
128 : }
129 : }
130 : });
131 :
132 : //
133 : // ELASTIC DRIVING FORCE
134 : //
135 30662 : if (pf.elastic_df)
136 : {
137 0 : amrex::Array4<const Set::Matrix> const& sigma = (*this->stress_mf[lev]).array(mfi);
138 :
139 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
140 : {
141 0 : for (int m = 0; m < number_of_grains; m++)
142 : {
143 0 : Set::Matrix F0avg = Set::Matrix::Zero();
144 :
145 0 : for (int n = 0; n < number_of_grains; n++)
146 : {
147 0 : F0avg += eta(i, j, k, n) * mechanics.model[n].F0;
148 : }
149 :
150 0 : 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 0 : Set::Matrix dF0deta = Set::Matrix::Zero();
154 :
155 0 : for (int n = 0; n < number_of_grains; n++)
156 : {
157 0 : if (n == m) continue;
158 0 : Set::Scalar normsq = eta(i, j, k, m) * eta(i, j, k, m) + eta(i, j, k, n) * eta(i, j, k, n);
159 0 : 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 0 : Set::Scalar tmpdf = (dF0deta.transpose() * sig).trace();
164 :
165 0 : if (pf.threshold.mechanics)
166 0 : driving_force_threshold(i, j, k, m) -= pf.elastic_mult * tmpdf;
167 : else
168 0 : 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 30662 : if (!anisotropic_kinetics.on || time < anisotropic_kinetics.tstart)
180 : {
181 30047762 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
182 : {
183 10005700 : for (int m = 0; m < number_of_grains; m++)
184 : {
185 8007300 : if (pf.threshold.on)
186 : {
187 0 : if (driving_force_threshold(i, j, k, m) > pf.threshold.value)
188 : {
189 0 : if (pf.threshold.type == ThresholdType::Continuous)
190 0 : eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m) - pf.threshold.value);
191 0 : else if (pf.threshold.type == ThresholdType::Chop)
192 0 : eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m));
193 : }
194 0 : else if (driving_force_threshold(i, j, k, m) < -pf.threshold.value)
195 : {
196 0 : if (pf.threshold.type == ThresholdType::Continuous)
197 0 : eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m) + pf.threshold.value);
198 0 : else if (pf.threshold.type == ThresholdType::Chop)
199 0 : eta(i, j, k, m) -= pf.L * dt * (driving_force_threshold(i, j, k, m));
200 : }
201 : }
202 :
203 24021900 : 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 324 : if (anisotropic_kinetics.on && time >= anisotropic_kinetics.tstart)
215 : {
216 0 : for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
217 : {
218 0 : amrex::Box bx = mfi.tilebox();
219 0 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
220 0 : amrex::Array4<Set::Scalar> const& L = (*anisotropic_kinetics.L_mf[lev]).array(mfi);
221 0 : amrex::Array4<Set::Scalar> const& threshold = (*anisotropic_kinetics.threshold_mf[lev]).array(mfi);
222 :
223 0 : for (int m = 0; m < number_of_grains; m++)
224 : {
225 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
226 : {
227 0 : Set::Vector Deta = Numeric::Gradient(eta, i, j, k, m, DX);
228 0 : Set::Scalar theta = atan2(Deta(1), Deta(0));
229 0 : L(i, j, k, m) = (4. / 3.) * anisotropic_kinetics.mobility(theta) / pf.l_gb;
230 0 : threshold(i, j, k, m) = anisotropic_kinetics.threshold(theta);
231 : });
232 : }
233 : }
234 :
235 0 : for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
236 : {
237 0 : amrex::Box bx = mfi.tilebox();
238 0 : Set::Patch<Set::Scalar> eta = eta_mf.Patch(lev,mfi);
239 0 : Set::Patch<const Set::Scalar> driving_force = driving_force_mf.Patch(lev,mfi);
240 0 : Set::Patch<const Set::Scalar> driving_force_threshold = driving_force_threshold_mf.Patch(lev,mfi);
241 0 : Set::Patch<const Set::Scalar> L = anisotropic_kinetics.L_mf.Patch(lev,mfi);
242 0 : Set::Patch<const Set::Scalar> threshold = anisotropic_kinetics.threshold_mf.Patch(lev,mfi);
243 :
244 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
245 : {
246 0 : for (int m = 0; m < number_of_grains; m++)
247 : {
248 0 : if (pf.threshold.on)
249 : {
250 0 : if (driving_force_threshold(i, j, k, m) > threshold(i,j,k,m))
251 : {
252 0 : if (pf.threshold.type == ThresholdType::Continuous)
253 0 : eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m) - threshold(i,j,k,m));
254 0 : else if (pf.threshold.type == ThresholdType::Chop)
255 0 : eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m));
256 : }
257 0 : else if (driving_force_threshold(i, j, k, m) < -pf.threshold.value)
258 : {
259 0 : if (pf.threshold.type == ThresholdType::Continuous)
260 0 : eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m) + threshold(i,j,k,m));
261 0 : else if (pf.threshold.type == ThresholdType::Chop)
262 0 : eta(i, j, k, m) -= L(i,j,k,m) * dt * (driving_force_threshold(i, j, k, m));
263 : }
264 : }
265 :
266 0 : eta(i, j, k, m) -= L(i,j,k,m) * dt * driving_force(i, j, k, m);
267 : }
268 : });
269 : }
270 : }
271 324 : }
272 :
273 : template<class model_type>
274 23 : void PhaseFieldMicrostructure<model_type>::Initialize(int lev)
275 : {
276 : BL_PROFILE("PhaseFieldMicrostructure::Initialize");
277 23 : Base::Mechanics<model_type>::Initialize(lev);
278 23 : ic->Initialize(lev, eta_mf);
279 23 : }
280 :
281 : template<class model_type>
282 50 : void PhaseFieldMicrostructure<model_type>::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
283 : {
284 : BL_PROFILE("PhaseFieldMicrostructure::TagCellsForRefinement");
285 50 : Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
286 50 : const amrex::Real* DX = this->geom[lev].CellSize();
287 50 : const Set::Vector dx(DX);
288 50 : const Set::Scalar dxnorm = dx.lpNorm<2>();
289 :
290 3436 : for (amrex::MFIter mfi(*eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
291 : {
292 3386 : const amrex::Box& bx = mfi.tilebox();
293 3386 : amrex::Array4<const amrex::Real> const& etanew = (*eta_mf[lev]).array(mfi);
294 3386 : amrex::Array4<char> const& tags = a_tags.array(mfi);
295 :
296 18238 : for (int n = 0; n < number_of_grains; n++)
297 1199942 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
298 1185090 : Set::Vector grad = Numeric::Gradient(etanew, i, j, k, n, DX);
299 :
300 1185090 : if (dxnorm * grad.lpNorm<2>() > ref_threshold)
301 115862 : tags(i, j, k) = amrex::TagBox::SET;
302 : });
303 : }
304 50 : }
305 :
306 : template<class model_type>
307 44 : void PhaseFieldMicrostructure<model_type>::TimeStepComplete(Set::Scalar /*time*/, int /*iter*/)
308 : {
309 : // TODO: remove this function, it is no longer needed.
310 44 : }
311 :
312 : template<class model_type>
313 0 : void PhaseFieldMicrostructure<model_type>::UpdateModel(int a_step, Set::Scalar /*a_time*/)
314 : {
315 : BL_PROFILE("PhaseFieldMicrostructure::UpdateModel");
316 0 : if (a_step % this->m_interval) return;
317 :
318 0 : for (int lev = 0; lev <= this->finest_level; ++lev)
319 : {
320 0 : amrex::Box domain = this->geom[lev].Domain();
321 0 : domain.convert(amrex::IntVect::TheNodeVector());
322 :
323 0 : eta_mf[lev]->FillBoundary();
324 :
325 0 : for (MFIter mfi(*this->model_mf[lev], false); mfi.isValid(); ++mfi)
326 : {
327 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
328 :
329 0 : amrex::Array4<model_type> const& model = this->model_mf[lev]->array(mfi);
330 0 : amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
331 :
332 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
333 : {
334 0 : std::vector<Set::Scalar> etas(number_of_grains);
335 0 : for (int n = 0; n < number_of_grains; n++)
336 0 : etas[n] = Numeric::Interpolate::CellToNodeAverage(eta, i, j, k, n);
337 0 : model(i, j, k) = model_type::Combine(mechanics.model, etas, 1);
338 : });
339 : }
340 :
341 0 : Util::RealFillBoundary(*this->model_mf[lev], this->geom[lev]);
342 : }
343 :
344 : }
345 :
346 : template<class model_type>
347 44 : void PhaseFieldMicrostructure<model_type>::TimeStepBegin(Set::Scalar time, int iter)
348 : {
349 : BL_PROFILE("PhaseFieldMicrostructure::TimeStepBegin");
350 44 : Base::Mechanics<model_type>::TimeStepBegin(time, iter);
351 :
352 44 : if (anisotropy.on && time >= anisotropy.tstart)
353 : {
354 40 : this->SetTimestep(anisotropy.timestep);
355 40 : if (anisotropy.elastic_int > 0)
356 0 : if (iter % anisotropy.elastic_int) return;
357 : }
358 : }
359 :
360 : template<class model_type>
361 11898 : 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 11898 : Base::Mechanics<model_type>::Integrate(amrlev, time, step, mfi, box);
366 :
367 11898 : Model::Interface::GB::SH gbmodel(0.0, 0.0, anisotropy.sigma0, anisotropy.sigma1);
368 11898 : const amrex::Real* DX = this->geom[amrlev].CellSize();
369 11898 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
370 :
371 11898 : amrex::Array4<amrex::Real> const& eta = (*eta_mf[amrlev]).array(mfi);
372 574826 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
373 : #if AMREX_SPACEDIM == 2
374 562928 : auto sten = Numeric::GetStencil(i, j, k, box);
375 : #endif
376 :
377 562928 : volume += eta(i, j, k, 0) * dv;
378 :
379 562928 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
380 562928 : Set::Scalar normgrad = grad.lpNorm<2>();
381 :
382 562928 : if (normgrad > 1E-8)
383 : {
384 169690 : Set::Vector normal = grad / normgrad;
385 :
386 169690 : Set::Scalar da = normgrad * dv;
387 169690 : area += da;
388 :
389 169690 : if (!anisotropy.on || time < anisotropy.tstart)
390 : {
391 8562 : gbenergy += pf.sigma0 * da;
392 8562 : Set::Scalar k = 0.75 * pf.sigma0 * pf.l_gb;
393 8562 : realgbenergy += 0.5 * k * normgrad * normgrad * dv;
394 8562 : regenergy = 0.0;
395 : }
396 : else
397 : {
398 : #if AMREX_SPACEDIM == 2
399 161128 : Set::Scalar theta = atan2(grad(1), grad(0));
400 161128 : Set::Scalar sigma = boundary->W(theta);
401 161128 : gbenergy += sigma * da;
402 :
403 161128 : Set::Scalar k = 0.75 * sigma * pf.l_gb;
404 161128 : realgbenergy += 0.5 * k * normgrad * normgrad * dv;
405 :
406 161128 : Set::Matrix DDeta = Numeric::Hessian(eta, i, j, k, 0, DX, sten);
407 161128 : Set::Vector tangent(normal[1], -normal[0]);
408 161128 : Set::Scalar k2 = (DDeta * tangent).dot(tangent);
409 161128 : regenergy += 0.5 * anisotropy.beta * k2 * k2;
410 : #elif AMREX_SPACEDIM == 3
411 0 : gbenergy += gbmodel.W(normal) * da;
412 : #endif
413 : }
414 : }
415 : });
416 11898 : }
417 :
418 : template class PhaseFieldMicrostructure<Model::Solid::Affine::Cubic>;
419 : template class PhaseFieldMicrostructure<Model::Solid::Affine::Hexagonal>;
420 :
421 :
422 : } // namespace Integrator
|
