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