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 9 : 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 :
58 : // whether to manually specify disconnection nucleation points
59 0 : pp.query("disconnection.fixed.on",value.fixed.on);
60 0 : if (value.fixed.on)
61 : {
62 : // array of x locations
63 0 : pp.queryarr("fixed.sitex",value.fixed.sitex);
64 : // array of y locations
65 0 : pp.queryarr("fixed.sitey",value.fixed.sitey);
66 : // array of order parameter number
67 0 : pp.queryarr("fixed.phases",value.fixed.phases);
68 : // time to appear
69 0 : pp.queryarr("fixed.time",value.fixed.time);
70 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.sitey.size()));
71 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.phases.size()));
72 0 : Util::Assert(INFO,TEST(value.fixed.sitex.size() == value.fixed.time.size()));
73 0 : value.fixed.done.resize(value.fixed.sitex.size(),false);
74 : }
75 : else
76 : {
77 0 : value.unif_dist = std::uniform_real_distribution<double>(0.0,1.0);
78 0 : value.int_dist = std::uniform_int_distribution<int>(0,1);
79 0 : value.rand_num_gen.seed(amrex::ParallelDescriptor::MyProc());
80 : }
81 :
82 : // verbosity
83 0 : pp_query_default("verbose",value.verbose,false);
84 0 : }
85 :
86 :
87 : /// This operates on an entire field, and manages all of the MPI
88 : /// communication necessary for consistent nucleation.
89 : void
90 0 : Nucleate( Set::Field<Set::Scalar> &eta_mf,
91 : std::vector<amrex::Geometry> &geom,
92 : Set::Scalar timestep,
93 : Set::Scalar time,
94 : int iter
95 : )
96 : {
97 0 : Util::Assert(INFO,TEST(eta_mf[0]->nComp() == 2), "This only works for 2 component phase fields");
98 :
99 0 : if (time < tstart) // wait until it's time to go
100 0 : return;
101 0 : if (iter % interval) // skip every [interval] timesteps
102 0 : return;
103 :
104 0 : sitex.clear();
105 0 : sitey.clear();
106 0 : phases.clear();
107 :
108 0 : int max_lev = eta_mf.finest_level;
109 :
110 0 : const Set::Scalar *DX = geom[max_lev].CellSize();
111 0 : Set::Scalar exponent = DX[0] * DX[0] * (timestep * interval) / tau_vol;
112 :
113 : // Determine the nucleation sites in the finest grid only
114 0 : for (amrex::MFIter mfi(*eta_mf[max_lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
115 : {
116 0 : const amrex::Box &bx = mfi.tilebox();
117 0 : amrex::Array4<amrex::Real> const &eta = (*eta_mf[max_lev]).array(mfi);
118 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
119 : {
120 0 : Set::Scalar E0 = 2.0*nucleation_energy;
121 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);
122 0 : Set::Scalar p = std::exp(-E0/(K_b*temp));
123 0 : Set::Scalar P = 1.0 - std::pow(1.0 - p,exponent);
124 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;
125 0 : Set::Scalar q = 0.0;
126 0 : q = unif_dist(rand_num_gen);
127 0 : if (q < P)
128 : {
129 0 : sitex.push_back(geom[max_lev].ProbLo()[0] + ((amrex::Real)(i)) * DX[0]);
130 0 : sitey.push_back(geom[max_lev].ProbLo()[1] + ((amrex::Real)(j)) * DX[1]);
131 0 : int phase = int_dist(rand_num_gen);
132 0 : phases.push_back(phase);
133 0 : } });
134 0 : }
135 :
136 : // Sync up all the nucleation sites among processors
137 0 : Util::MPI::Allgather(sitex);
138 0 : Util::MPI::Allgather(sitey);
139 0 : Util::MPI::Allgather(phases);
140 :
141 0 : if (verbose)
142 : {
143 0 : Util::Message(INFO, "Nucleating ", phases.size(), " disconnections");
144 : }
145 :
146 0 : if (sitex.size() > 0)
147 : {
148 : // Now that we all know the nucleation locations, perform the nucleation
149 0 : for (int lev = 0; lev <= max_lev; lev++)
150 : {
151 0 : amrex::Box domain = geom[lev].Domain();
152 0 : domain.convert(amrex::IntVect::TheNodeVector());
153 0 : const amrex::Real *DX = geom[lev].CellSize();
154 0 : for (amrex::MFIter mfi(*eta_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
155 : {
156 0 : const amrex::Box bx = mfi.grownnodaltilebox() & domain;
157 0 : amrex::Array4<Set::Scalar> const &eta = (*eta_mf[lev]).array(mfi);
158 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
159 : {
160 0 : Set::Vector x;
161 0 : AMREX_D_TERM(
162 : x(0) = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * DX[0];,
163 : x(1) = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * DX[1];,
164 : x(2) = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * DX[2];);
165 0 : for (unsigned int m = 0; m < phases.size(); m++)
166 : {
167 0 : amrex::Real r_squared = 0;
168 0 : Set::Vector nucleation_site(AMREX_D_DECL(sitex[m], sitey[m], 0.0));
169 0 : for (int n = 0; n < AMREX_SPACEDIM; n++)
170 : {
171 0 : amrex::Real dist = nucleation_site(n) - x(n);
172 0 : r_squared += dist * dist;
173 : }
174 0 : amrex::Real bump = exp(-r_squared / box_size);
175 0 : eta(i, j, k, phases[m]) = bump * (1 - eta(i, j, k, phases[m])) + eta(i, j, k, phases[m]);
176 0 : eta(i, j, k, 1 - phases[m]) = (1. - bump) * eta(i, j, k, 1 - phases[m]);
177 : }
178 0 : });
179 0 : }
180 : }
181 : }
182 : }
183 :
184 :
185 :
186 : private:
187 :
188 : bool verbose = false;
189 :
190 : Set::Scalar tstart = NAN;
191 : Set::Scalar nucleation_energy = NAN;
192 : Set::Scalar tau_vol = NAN;
193 : Set::Scalar temp = NAN;
194 : Set::Scalar box_size = NAN;
195 : Set::Scalar epsilon = NAN;
196 :
197 : int interval = -1;
198 :
199 : std::uniform_real_distribution<double> unif_dist; /// random number distribution for spatial location
200 : std::uniform_int_distribution<int> int_dist; /// random number generator for phase
201 : std::default_random_engine rand_num_gen; /// generator object
202 :
203 : struct {
204 : int on = 0;
205 : std::vector<Set::Scalar> sitex;
206 : std::vector<Set::Scalar> sitey;
207 : std::vector<int> phases;
208 : std::vector<Set::Scalar> time;
209 : std::vector<bool> done;
210 : } fixed;
211 :
212 :
213 : std::vector<Set::Scalar> sitex; /// list of nucleation site x coordinates
214 : std::vector<Set::Scalar> sitey; /// list of nucleation stie y coordinates
215 : std::vector<int> phases; /// list of nucleation site phases (up or down)
216 :
217 : const Set::Scalar K_b = 8.617333262145e-5; // eV/K
218 : };
219 : }
220 : }
221 :
222 : #endif
|
