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