Line data Source code
1 : ///
2 : /// \file PhaseFieldMicrostructure.H
3 : ///
4 : #ifndef INTEGRATOR_PHASEFIELDMICROSTRUCTURE_H
5 : #define INTEGRATOR_PHASEFIELDMICROSTRUCTURE_H
6 :
7 : #include <iostream>
8 : #include <fstream>
9 : #include <iomanip>
10 : #include <random>
11 : #include <chrono>
12 :
13 : #include "AMReX.H"
14 : #include "AMReX_ParmParse.H"
15 : #include "AMReX_ParallelDescriptor.H"
16 : #include <AMReX_MLMG.H>
17 :
18 : #include "Integrator/Integrator.H"
19 : #include "Integrator/Mechanics.H"
20 :
21 : #include "BC/BC.H"
22 : #include "BC/Constant.H"
23 : #include "BC/Operator/Elastic/Constant.H"
24 : #include "IC/Constant.H"
25 : #include "IC/TabulatedInterface.H"
26 : #include "IC/PerturbedInterface.H"
27 : #include "IC/Voronoi.H"
28 : #include "IC/Sphere.H"
29 : #include "IC/Expression.H"
30 : #include "IC/Ellipse.H"
31 : #include "IC/Random.H"
32 :
33 : #include "Model/Interface/GB/GB.H"
34 : #include "Model/Interface/GB/Sin.H"
35 : #include "Model/Interface/GB/AbsSin.H"
36 : #include "Model/Interface/GB/Read.H"
37 : #include "Model/Interface/GB/SH.H"
38 :
39 : #include "Model/Solid/Linear/Cubic.H"
40 : #include "Model/Solid/Affine/Cubic.H"
41 :
42 : #include "Model/Defect/Disconnection.H"
43 :
44 : #include "Operator/Elastic.H"
45 :
46 : namespace Integrator
47 : {
48 : //using model_type = Model::Solid::Affine::Cubic;
49 :
50 : enum RegularizationType{
51 : Wilmore = 1,
52 : K12 = 2};
53 :
54 : ///
55 : /// \class PhaseFieldMicrostructure
56 : /// \brief Microstructure evolution with grain boundary anisotropy
57 : ///
58 : /// Solve the Allen-Cahn evolution equation for microstructure with parameters \f$\eta_1\ldots\eta_n\f$,
59 : /// where n corresponds to the number of grains.
60 : ///
61 : template<class model_type>
62 : class PhaseFieldMicrostructure : public Base::Mechanics<model_type>
63 : {
64 : public:
65 0 : PhaseFieldMicrostructure() : Integrator() {};
66 9 : PhaseFieldMicrostructure(IO::ParmParse &pp) : Integrator()
67 36 : {pp_queryclass(*this);}
68 6 : virtual ~PhaseFieldMicrostructure()
69 : {
70 2 : delete boundary;
71 3 : delete ic;
72 3 : delete mybc;
73 9 : }
74 9 : static void Parse(PhaseFieldMicrostructure &value, IO::ParmParse &pp)
75 : {
76 : BL_PROFILE("PhaseFieldMicrostructure::Parse");
77 :
78 : // Number of grain fields (may be more if using different IC)
79 54 : pp_query_default("pf.number_of_grains", value.number_of_grains,2);
80 54 : pp_query_required("pf.M", value.pf.M); // Mobility
81 : //pp_query_required("pf.mu", value.pf.mu); // Phase field :math:`\mu`
82 54 : pp_query_required("pf.gamma", value.pf.gamma); // Phase field :math:`\gamma`
83 54 : pp_query_required("pf.sigma0", value.pf.sigma0); // Initial GB energy if not using anisotropy
84 54 : pp_query_required("pf.l_gb", value.pf.l_gb); // Mobility
85 54 : pp_query_default("pf.elastic_df",value.pf.elastic_df,false); // Determine whether to use elastic driving force
86 54 : pp_query_default("pf.elastic_mult",value.pf.elastic_mult,1.0); // Multiplier of elastic energy
87 :
88 63 : pp_query_default("pf.threshold.value",value.pf.threshold.value,0.0); // Value used for thresholding kinetic relation
89 54 : pp_query_default("pf.threshold.chempot",value.pf.threshold.chempot,false); // Whether to include chemical potential in threshold
90 54 : pp_query_default("pf.threshold.boundary",value.pf.threshold.boundary,false); // Whether to include boundary energy in threshold
91 54 : pp_query_default("pf.threshold.corner",value.pf.threshold.corner,false); // Whether to include corner regularization in threshold
92 54 : pp_query_default("pf.threshold.lagrange",value.pf.threshold.lagrange,false); // Whether to include lagrange multiplier in threshold
93 45 : pp_query_default("pf.threshold.mechanics",value.pf.threshold.mechanics,false);// Whether to include mechanical driving force in threshold
94 : {
95 9 : value.pf.threshold.on =
96 9 : value.pf.threshold.chempot || value.pf.threshold.boundary ||
97 26 : value.pf.threshold.corner || value.pf.threshold.lagrange ||
98 8 : value.pf.threshold.mechanics;
99 :
100 18 : std::string type = "continuous";
101 63 : pp_query_validate("pf.threshold.type",type,{"continuous","chop"}); // Type of thresholding to use
102 9 : if (type == "continuous") value.pf.threshold.type = ThresholdType::Continuous;
103 1 : else if (type == "chop") value.pf.threshold.type = ThresholdType::Chop;
104 9 : }
105 :
106 9 : value.pf.L = (4./3.)*value.pf.M / value.pf.l_gb;
107 :
108 54 : pp_query_required("amr.max_level", value.max_level); // Maximum AMR level
109 54 : pp_query_default("amr.ref_threshold", value.ref_threshold, 0.1); // Phase field refinement threshold
110 :
111 : // Elasticity
112 45 : pp_query_default("mechanics.tstart",value.mechanics.tstart, 0.0);
113 9 : value.mechanics.model.clear();
114 9 : value.mechanics.model.resize(value.number_of_grains,model_type::Zero());
115 59 : for (int i=0; i < value.number_of_grains; i++)
116 : {
117 : // By default, read in the model specified by "mechanics.model"
118 100 : if (pp.getEntries("mechanics.model").size())
119 120 : pp_queryclass("mechanics.model", value.mechanics.model[i]);
120 :
121 : // If individual models are specified, read those in and overwrite
122 50 : std::string name = "mechanics.model" + std::to_string(i+1);
123 50 : if (pp.getEntries(name).size())
124 12 : pp_queryclass(std::string(name.data()), value.mechanics.model[i]);
125 : }
126 : // Mixing order
127 72 : pp_query_validate("mechanics.mix_order",value.mechanics.model_mix_order,{1,2});
128 : // Force Neumann BCs on the model
129 45 : pp_query_default("mechanics.model_neuman_boundary",value.mechanics.model_neumann_boundary,false);
130 :
131 9 : value.m_type = Base::Mechanics<model_type>::Disable; // Turn mechanics off by default
132 45 : pp.queryclass<Base::Mechanics<model_type>>("mechanics",value);
133 9 : if (value.m_type == Base::Mechanics<model_type>::Type::Static)
134 3 : value.number_of_ghost_cells = std::max(value.number_of_ghost_cells,3);
135 :
136 :
137 : // Lagrange multiplier method for enforcing volumes
138 45 : pp_query_default("lagrange.on", value.lagrange.on,false);
139 9 : if (value.lagrange.on)
140 : {
141 12 : pp_query_required("lagrange.lambda", value.lagrange.lambda); // Lagrange multiplier value
142 12 : pp_query_required("lagrange.vol0", value.lagrange.vol0); // Prescribed volume
143 10 : pp_query_default("lagrange.tstart", value.lagrange.tstart,0.0); // Lagrange multipler start time
144 2 : value.SetThermoInt(1);
145 : }
146 :
147 45 : pp_query_default("sdf.on",value.sdf.on,false); // synthetic driving force (SDF)
148 9 : if(value.sdf.on)
149 : {
150 0 : std::vector<std::string> vals;
151 0 : pp_queryarr("sdf.val",vals); // value of SDF for each grain
152 0 : int nvals = static_cast<int>(vals.size());
153 0 : if (nvals == 1)
154 0 : for (int i = 0; i < value.number_of_grains; i++)
155 0 : value.sdf.val.push_back(Numeric::Interpolator::Linear<Set::Scalar>(vals[0]));
156 0 : else if (nvals == value.number_of_grains)
157 0 : for (int i = 0; i < value.number_of_grains; i++)
158 0 : value.sdf.val.push_back(Numeric::Interpolator::Linear<Set::Scalar>(vals[i]));
159 : else
160 0 : Util::Abort(INFO,"sdf.val received ", vals.size(), " but requires 1 or ", value.number_of_grains);
161 :
162 0 : pp_query_default("sdf.tstart",value.sdf.tstart,0.0); // time to begin applying SDF
163 0 : }
164 :
165 : // Anisotropic grain boundary energy parameters
166 :
167 45 : pp_query_default("anisotropy.on", value.anisotropy.on,false); // Turn on
168 9 : if (value.anisotropy.on)
169 : {
170 : // Regularization para m
171 30 : pp_query_required("anisotropy.beta", value.anisotropy.beta);
172 : // Time to turn on anisotropy
173 25 : pp_query_required("anisotropy.tstart", value.anisotropy.tstart);
174 5 : value.anisotropy.timestep = value.timestep;
175 : // Modify timestep when turned on
176 25 : pp_query_required("anisotropy.timestep", value.anisotropy.timestep);
177 5 : value.anisotropy.plot_int = value.plot_int;
178 : // Modify plot_int when turned on
179 25 : pp_query_default("anisotropy.plot_int", value.anisotropy.plot_int, -1);
180 5 : value.anisotropy.plot_dt = value.plot_dt;
181 : // Modify plot_dt when turned on
182 30 : pp_query_default("anisotropy.plot_dt", value.anisotropy.plot_dt, -1.0);
183 : // Modify thermo int when turned on
184 30 : pp_query_default("anisotropy.thermo_int", value.anisotropy.thermo_int, -1);
185 : // Modify thermo plot int when turned on
186 30 : pp_query_default("anisotropy.thermo_plot_int", value.anisotropy.thermo_plot_int, -1);
187 : // Frequency of elastic calculation
188 25 : pp_query_default("anisotropy.elastic_int",value.anisotropy.elastic_int, -1);
189 5 : if (value.anisotropy.on)
190 5 : value.number_of_ghost_cells = std::max(value.number_of_ghost_cells,2);
191 :
192 : // Determine the kind of regularization to use
193 5 : std::map<std::string, RegularizationType> regularization_type;
194 10 : regularization_type["wilmore"] = RegularizationType::Wilmore;
195 5 : regularization_type["k12"] = RegularizationType::K12;
196 :
197 5 : std::string regularization_type_input;
198 : // Type of regularization to use
199 35 : pp_query_validate("anisotropy.regularization", regularization_type_input,{"k12","wilmore"});
200 5 : value.regularization = regularization_type[regularization_type_input];
201 :
202 35 : pp_forbid("anisotropy.gb_type"," --> anisotropy.type");
203 : // Type of GB to use
204 10 : pp.select<Model::Interface::GB::AbsSin, Model::Interface::GB::Sin, Model::Interface::GB::Read, Model::Interface::GB::SH>("anisotropy",value.boundary);
205 5 : }
206 : // Thermal fluctuations
207 9 : pp.query("fluctuation.on",value.fluctuation.on);
208 9 : if (value.fluctuation.on)
209 : {
210 0 : pp.query("fluctuation.amp",value.fluctuation.amp); // fluctuation amplitude
211 0 : pp.query("fluctuation.sd",value.fluctuation.sd); // fluctuation stadard deviation
212 0 : pp.query("fluctuation.tstart", value.fluctuation.tstart); // time to start applying fluctuation
213 0 : value.fluctuation.norm_dist = std::normal_distribution<double>(0.0,value.fluctuation.sd);
214 : }
215 :
216 : // Disconnection generation
217 9 : pp.query("disconnection.on",value.disconnection.on);
218 9 : if (value.disconnection.on)
219 : {
220 : // Read in nucleation parameters from disconnection class
221 0 : pp.queryclass<Model::Defect::Disconnection>("disconnection",value.disconnection.model);
222 : }
223 :
224 : // Shear coupling matrices
225 45 : pp_query_default("shearcouple.on",value.shearcouple.on,false);
226 9 : if (value.shearcouple.on)
227 : {
228 0 : Util::AssertException(INFO,TEST(value.m_time_evolving == true), " mechanics.time_evolving must be true when using shearcouple");
229 0 : Util::AssertException(INFO,TEST(value.mechanics.model_mix_order == 2), " mechanics.model.mix_order must be 2 when using shearcouple");
230 :
231 0 : value.shearcouple.Fgb.resize(value.number_of_grains * value.number_of_grains, Set::Matrix::Zero());
232 0 : for (int i = 0 ; i < value.number_of_grains ; i++)
233 0 : for (int j = 0 ; j < value.number_of_grains ; j++)
234 : {
235 0 : std::string name = "shearcouple.Fgb."+std::to_string(i)+"."+std::to_string(j);
236 0 : std::string namerev = "shearcouple.Fgb."+std::to_string(j)+"."+std::to_string(i);
237 0 : if ( i==j && pp.contains(name.data()))
238 0 : Util::Abort(INFO,"Cannot specify self FGB ", name);
239 0 : if (pp.contains(name.data()) && pp.contains(namerev.data()))
240 0 : Util::Abort(INFO,"Cannot specify both ",name," and ",namerev);
241 0 : if (pp.contains(name.data()))
242 : {
243 0 : pp.queryarr(name.data(),value.shearcouple.Fgb[i*value.number_of_grains + j]);
244 0 : value.shearcouple.Fgb[j*value.number_of_grains + i] = - value.shearcouple.Fgb[i*value.number_of_grains + j];
245 : }
246 : }
247 : }
248 :
249 : // Boundary condition for eta
250 18 : pp.select<BC::Constant>("bc.eta",value.mybc,value.number_of_grains);
251 :
252 : // Initial condition for the order parameter eta
253 27 : pp.select<IC::Constant,IC::PerturbedInterface,IC::Voronoi,IC::Expression,IC::Sphere,IC::Ellipse,IC::Random>("ic",value.ic,value.geom);
254 :
255 : // Anisotropic mobility
256 45 : pp_query_default("anisotropic_kinetics.on",value.anisotropic_kinetics.on, 0);
257 9 : if (value.anisotropic_kinetics.on)
258 : {
259 : // simulation time when anisotropic kinetics is activated
260 5 : pp_query_default("anisotropic_kinetics.tstart",value.anisotropic_kinetics.tstart, 0.0);
261 1 : std::string mobility_filename, threshold_filename;
262 : // file containing anisotropic mobility data
263 5 : pp_query_file("anisotropic_kinetics.mobility",mobility_filename);
264 1 : value.anisotropic_kinetics.mobility = Numeric::Interpolator::Linear<Set::Scalar>::Read(mobility_filename);
265 : // file containing anisotropic mobility data
266 5 : pp_query_file("anisotropic_kinetics.threshold",threshold_filename);
267 1 : value.anisotropic_kinetics.threshold = Numeric::Interpolator::Linear<Set::Scalar>::Read(threshold_filename);
268 3 : value.RegisterNewFab(value.anisotropic_kinetics.L_mf, value.mybc, value.number_of_grains, 0, "mobility",true);
269 3 : value.RegisterNewFab(value.anisotropic_kinetics.threshold_mf, value.mybc, value.number_of_grains, 0, "theshold",true);
270 1 : }
271 :
272 :
273 :
274 :
275 27 : value.RegisterNewFab(value.eta_mf, value.mybc, value.number_of_grains, value.number_of_ghost_cells, "Eta",true);
276 27 : value.RegisterNewFab(value.eta_old_mf, value.mybc, value.number_of_grains, value.number_of_ghost_cells, "EtaOld",false);
277 : //value.RegisterNewFab(value.driving_force_mf, value.mybc, value.number_of_grains, value.number_of_ghost_cells, "DrivingForce",false);
278 : // if (value.pf.threshold.on)
279 : // value.RegisterNewFab(value.driving_force_threshold_mf, value.mybc, value.number_of_grains, value.number_of_ghost_cells, "DrivingForceThreshold",false);
280 : // if (value.disconnection.on)
281 : // value.RegisterNewFab(value.disc_mf, new BC::Nothing(), 1, value.number_of_ghost_cells, "disc",true); // see box
282 :
283 18 : value.RegisterIntegratedVariable(&value.volume, "volume");
284 18 : value.RegisterIntegratedVariable(&value.area, "area");
285 18 : value.RegisterIntegratedVariable(&value.gbenergy, "gbenergy");
286 18 : value.RegisterIntegratedVariable(&value.realgbenergy, "realgbenergy");
287 18 : value.RegisterIntegratedVariable(&value.regenergy, "regenergy");
288 :
289 9 : }
290 :
291 :
292 : protected:
293 :
294 : /// \fn Advance
295 : /// \brief Evolve phase field in time
296 : void Advance (int lev, Real time, Real dt) override;
297 : void Initialize (int lev) override;
298 :
299 : void TagCellsForRefinement (int lev, amrex::TagBoxArray& tags, amrex::Real time, int ngrow) override;
300 :
301 : virtual void TimeStepBegin(amrex::Real time, int iter) override;
302 : virtual void TimeStepComplete(amrex::Real time, int iter) override;
303 : void Integrate(int amrlev, Set::Scalar time, int step,
304 : const amrex::MFIter &mfi, const amrex::Box &box) override;
305 :
306 : virtual void UpdateEigenstrain(int lev);
307 0 : virtual void UpdateEigenstrain()
308 : {
309 0 : for (int lev = 0; lev <= this->max_level; lev++)
310 0 : UpdateEigenstrain(lev);
311 0 : }
312 :
313 : virtual void UpdateModel(int /*a_step*/, Set::Scalar /*a_time*/) override;
314 :
315 : private:
316 :
317 : int number_of_grains = -1;
318 : int number_of_ghost_cells = 1;
319 : Set::Scalar ref_threshold = 0.1;
320 :
321 : // Cell fab
322 : Set::Field<Set::Scalar> eta_mf; // Multicomponent field variable storing \t$\eta_i\t$ for the __current__ timestep
323 : Set::Field<Set::Scalar> eta_old_mf;
324 : Set::Field<Set::Scalar> driving_force_mf;
325 : Set::Field<Set::Scalar> driving_force_threshold_mf;
326 : Set::Field<Set::Scalar> fluct_mf;
327 : //Set::Field<Set::Scalar> disc_mf; //see box
328 : // Node fab
329 : //Set::Field<Set::Scalar> elasticdf_mf;
330 : Set::Field<Set::Scalar> totaldf_mf;
331 :
332 : BC::BC<Set::Scalar> *mybc = nullptr;
333 :
334 : //amrex::Real M, mu, gamma, sigma0, l_gb, beta;
335 : RegularizationType regularization = RegularizationType::K12;
336 : enum ThresholdType {
337 : Continuous = 0, Chop = 1
338 : };
339 : struct {
340 : Set::Scalar M = NAN;
341 : Set::Scalar L = NAN;
342 : Set::Scalar mu = NAN;
343 : Set::Scalar gamma = NAN;
344 : Set::Scalar sigma0 = NAN;
345 : Set::Scalar l_gb = NAN;
346 : bool elastic_df = false;
347 : Set::Scalar elastic_mult = NAN;
348 : struct {
349 : bool on = false;
350 : bool chempot = false;
351 : bool boundary = false;
352 : bool corner = false;
353 : bool lagrange = false;
354 : bool mechanics = false;
355 : bool sdf = false;
356 : Set::Scalar value = NAN;
357 : ThresholdType type = ThresholdType::Continuous;
358 : } threshold;
359 : } pf;
360 :
361 : struct {
362 : int on = 0;
363 : Set::Scalar tstart = 0.0;
364 : Numeric::Interpolator::Linear<Set::Scalar> mobility;
365 : Numeric::Interpolator::Linear<Set::Scalar> threshold;
366 : Set::Field<Set::Scalar> L_mf;
367 : Set::Field<Set::Scalar> threshold_mf;
368 : } anisotropic_kinetics;
369 :
370 : struct {
371 : int on = 0;
372 : Set::Scalar beta;
373 : Set::Scalar timestep;
374 : Set::Scalar tstart;
375 : int plot_int = -1;
376 : Set::Scalar plot_dt = -1.0;
377 : int thermo_int = -1, thermo_plot_int = -1;
378 : Set::Scalar theta0,sigma0,sigma1;
379 : Set::Scalar phi0 = 0.0;
380 : int elastic_int = -1;
381 : } anisotropy;
382 :
383 : struct {
384 : bool on = 0;
385 : Set::Scalar tstart = NAN;
386 : Set::Scalar vol0 = NAN;
387 : Set::Scalar lambda = NAN;
388 : } lagrange;
389 :
390 : struct {
391 : int on = 0;
392 : std::vector<Numeric::Interpolator::Linear<Set::Scalar>> val;
393 : Set::Scalar tstart = 0.0;
394 : } sdf;
395 :
396 : struct {
397 : int on = 0;
398 : Set::Scalar amp = 0.0;
399 : Set::Scalar sd = 0.0;
400 : Set::Scalar tstart = 0.0;
401 : std::normal_distribution<double> norm_dist;
402 : std::default_random_engine rand_num_gen;
403 : } fluctuation;
404 :
405 : struct {
406 : int on = 0;
407 : Model::Defect::Disconnection model;
408 : } disconnection;
409 :
410 : struct {
411 : int on = 0;
412 : std::vector<Set::Matrix> Fgb;
413 : } shearcouple;
414 :
415 : std::string gb_type, filename;
416 :
417 : Model::Interface::GB::GB *boundary = nullptr;
418 :
419 : IC::IC<Set::Scalar> *ic = nullptr;
420 :
421 : Set::Scalar volume = 5;
422 : Set::Scalar area = 0.0;
423 : Set::Scalar gbenergy = 0.0;
424 : Set::Scalar realgbenergy = 0.0;
425 : Set::Scalar regenergy = 0.0;
426 :
427 : struct
428 : {
429 : Set::Scalar tstart = 0.0;
430 : std::vector<model_type> model;
431 : int model_mix_order;
432 : bool model_neumann_boundary = false;
433 : } mechanics;
434 :
435 : using Base::Mechanics<model_type>::model_mf;
436 : using Base::Mechanics<model_type>::stress_mf;
437 :
438 : };
439 : }
440 : #endif
|
