Line data Source code
1 : //
2 : // This simulates the creation of "disconnection pairs" at a GB
3 : // by perturbing the order parameter :math:`\eta` with a gaussian.
4 : //
5 : // Disconnections can be nucleated randomly, or they can be added
6 : // in a controlled way through the "fixed" option.
7 : //
8 : // See the following reference for further details
9 : //
10 : // .. bibliography::
11 : // :list: none
12 : // :filter: False
13 : //
14 : // gokuli2021multiphase
15 : //
16 : //
17 :
18 : #ifndef MODEL_DEFECT_DISCONNECTION
19 : #define MODEL_DEFECT_DISCONNECTION
20 :
21 : #include <random>
22 :
23 : #include "AMReX_SPACE.H"
24 : #include "IO/ParmParse.H"
25 : #include "Set/Base.H"
26 : #include "Util/MPI.H"
27 :
28 : namespace Model
29 : {
30 : namespace Defect
31 : {
32 : class Disconnection
33 : {
34 : public:
35 :
36 3 : Disconnection(){};
37 3 : ~Disconnection() {};
38 :
39 : static void
40 0 : Parse(Disconnection &value, IO::ParmParse &pp)
41 : {
42 0 : Util::Assert(INFO,TEST(AMREX_SPACEDIM==2),"2D only");
43 : // time to start applying disconnections
44 0 : pp_query_default ("tstart", value.tstart,0.0);
45 : // nucleation energy
46 0 : pp_query_default ("nucleation_energy", value.nucleation_energy,0.0);
47 : // characteristic time
48 0 : pp_query_default("tau_vol", value.tau_vol, 1.0);
49 : // temperature
50 0 : pp_query_default("temp", value.temp, 0.0);
51 : // characteristic size
52 0 : pp_query_default("box_size", value.box_size, 0.0);
53 : // interval between generation events
54 0 : pp_query_required("interval", value.interval);
55 : // regularization epsilon
56 0 : pp_query_default ("epsilon",value.epsilon,1E-20);
57 : // whether to manually specify disconnection nucleation points
58 0 : pp.query("disconnection.fixed.on",value.fixed.on);
59 0 : if (value.fixed.on)
60 : {
61 : // array of x locations
62 0 : pp.queryarr("fixed.sitex",value.fixed.sitex);
63 : // array of y locations
64 0 : pp.queryarr("fixed.sitey",value.fixed.sitey);
65 : // array of order parameter number
66 0 : pp.queryarr("fixed.phases",value.fixed.phases);
67 : // time to appear
68 0 : pp.queryarr("fixed.time",value.fixed.time);
69 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.sitey.size()));
70 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.phases.size()));
71 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.time.size()));
72 0 : value.fixed.done.resize(value.fixed.sitex.size(),false);
73 : }
74 : else
75 : {
76 0 : value.unif_dist = std::uniform_real_distribution<double>(0.0,1.0);
77 0 : value.int_dist = std::uniform_int_distribution<int>(0,1);
78 0 : value.rand_num_gen.seed(amrex::ParallelDescriptor::MyProc());
79 : }
80 :
81 : // verbosity
82 0 : pp_query_default("verbose",value.verbose,false);
83 0 : }
84 :
85 :
86 : /// This operates on an entire field, and manages all of the MPI
87 : /// communication necessary for consistent nucleation.
88 : void
89 0 : Nucleate( Set::Field<Set::Scalar> &eta_mf,
90 : std::vector<amrex::Geometry> &geom,
91 : Set::Scalar timestep,
92 : Set::Scalar time,
93 : int iter
94 : )
95 : {
96 0 : Util::Assert(INFO,TEST(eta_mf[0]->nComp() == 2), "This only works for 2 component phase fields");
97 :
98 0 : if (time < tstart) // wait until it's time to go
99 0 : return;
100 0 : if (iter % interval) // skip every [interval] timesteps
101 0 : return;
102 :
103 0 : sitex.clear();
104 0 : sitey.clear();
105 0 : phases.clear();
106 :
107 0 : int max_lev = eta_mf.finest_level;
108 :
109 0 : const Set::Scalar *DX = geom[max_lev].CellSize();
110 0 : Set::Scalar exponent = DX[0] * DX[0] * (timestep * interval) / tau_vol;
111 :
112 : // Determine the nucleation sites in the finest grid only
113 0 : for (amrex::MFIter mfi(*eta_mf[max_lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
114 : {
115 0 : const amrex::Box &bx = mfi.tilebox();
116 0 : amrex::Array4<amrex::Real> const &eta = (*eta_mf[max_lev]).array(mfi);
117 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
118 : {
119 0 : Set::Scalar E0 = 2.0*nucleation_energy;
120 0 : E0 /= epsilon + 256.0*eta(i,j,k,0)*eta(i,j,k,0)*eta(i,j,k,0)*eta(i,j,k,0)*eta(i,j,k,1)*eta(i,j,k,1)*eta(i,j,k,1)*eta(i,j,k,1);
121 0 : Set::Scalar p = std::exp(-E0/(K_b*temp));
122 0 : Set::Scalar P = 1.0 - std::pow(1.0 - p,exponent);
123 0 : if (eta(i,j,k,0) < 0 || eta(i,j,k,0) > 1.0 || eta(i,j,k,1) < 0 || eta(i,j,k,1) > 1.0) P = 0.0;
124 0 : Set::Scalar q = 0.0;
125 0 : q = unif_dist(rand_num_gen);
126 0 : if (q < P)
127 : {
128 0 : sitex.push_back(geom[max_lev].ProbLo()[0] + ((amrex::Real)(i)) * DX[0]);
129 0 : sitey.push_back(geom[max_lev].ProbLo()[1] + ((amrex::Real)(j)) * DX[1]);
130 0 : int phase = int_dist(rand_num_gen);
131 0 : phases.push_back(phase);
132 0 : } });
133 0 : }
134 :
135 : // Sync up all the nucleation sites among processors
136 0 : Util::MPI::Allgather(sitex);
137 0 : Util::MPI::Allgather(sitey);
138 0 : Util::MPI::Allgather(phases);
139 :
140 0 : if (verbose)
141 : {
142 0 : Util::Message(INFO, "Nucleating ", phases.size(), " disconnections");
143 : }
144 :
145 0 : if (sitex.size() > 0)
146 : {
147 : // Now that we all know the nucleation locations, perform the nucleation
148 0 : for (int lev = 0; lev <= max_lev; lev++)
149 : {
150 0 : amrex::Box domain = geom[lev].Domain();
151 0 : domain.convert(amrex::IntVect::TheNodeVector());
152 0 : const amrex::Real *DX = geom[lev].CellSize();
153 0 : for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
154 : {
155 0 : const amrex::Box bx = mfi.grownnodaltilebox() & domain;
156 0 : amrex::Array4<Set::Scalar> const &eta = (*eta_mf[lev]).array(mfi);
157 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
158 : {
159 0 : Set::Vector x;
160 0 : AMREX_D_TERM(
161 : x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * DX[0];,
162 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * DX[1];,
163 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * DX[2];);
164 0 : for (unsigned int m = 0; m < phases.size(); m++)
165 : {
166 0 : amrex::Real r_squared = 0;
167 0 : Set::Vector nucleation_site(AMREX_D_DECL(sitex[m], sitey[m], 0.0));
168 0 : for (int n = 0; n < AMREX_SPACEDIM; n++)
169 : {
170 0 : amrex::Real dist = nucleation_site(n) - x(n);
171 0 : r_squared += dist * dist;
172 : }
173 0 : amrex::Real bump = exp(-r_squared / box_size);
174 0 : eta(i, j, k, phases[m]) = bump * (1 - eta(i, j, k, phases[m])) + eta(i, j, k, phases[m]);
175 0 : eta(i, j, k, 1 - phases[m]) = (1. - bump) * eta(i, j, k, 1 - phases[m]);
176 : }
177 0 : });
178 0 : }
179 : }
180 : }
181 : }
182 :
183 :
184 :
185 : private:
186 :
187 : bool verbose = false;
188 :
189 : Set::Scalar tstart = NAN;
190 : Set::Scalar nucleation_energy = NAN;
191 : Set::Scalar tau_vol = NAN;
192 : Set::Scalar temp = NAN;
193 : Set::Scalar box_size = NAN;
194 : Set::Scalar epsilon = NAN;
195 :
196 : int interval = -1;
197 :
198 : std::uniform_real_distribution<double> unif_dist; /// random number distribution for spatial location
199 : std::uniform_int_distribution<int> int_dist; /// random number generator for phase
200 : std::default_random_engine rand_num_gen; /// generator object
201 :
202 : struct {
203 : int on = 0;
204 : std::vector<Set::Scalar> sitex;
205 : std::vector<Set::Scalar> sitey;
206 : std::vector<int> phases;
207 : std::vector<Set::Scalar> time;
208 : std::vector<bool> done;
209 : } fixed;
210 :
211 :
212 : std::vector<Set::Scalar> sitex; /// list of nucleation site x coordinates
213 : std::vector<Set::Scalar> sitey; /// list of nucleation stie y coordinates
214 : std::vector<int> phases; /// list of nucleation site phases (up or down)
215 :
216 : const Set::Scalar K_b = 8.617333262145e-5; // eV/K
217 : };
218 : }
219 : }
220 :
221 : #endif
|
