Line data Source code
1 : //
2 : // This is a BASIC phase field method for dendritic growth as described in
3 : // "A Numerical Approach to Three-Dimensional Dendritic Solidification" by Ryo Kobayashi (1994)
4 : // Reference: https://doi.org/10.1080/10586458.1994.10504577
5 : //
6 :
7 : #ifndef INTEGRATOR_DENDRITE_H // Include guards
8 : #define INTEGRATOR_DENDRITE_H //
9 :
10 : #include <iostream>
11 : #include <fstream>
12 : #include <iomanip>
13 :
14 : #include "AMReX.H"
15 : #include "AMReX_ParallelDescriptor.H"
16 : #include "AMReX_ParmParse.H"
17 :
18 : #include "IO/ParmParse.H"
19 : #include "Integrator/Integrator.H"
20 : #include "BC/Constant.H"
21 : #include "IC/IC.H"
22 : #include "IC/Sphere.H"
23 : #include "IC/Constant.H"
24 : #include "IC/Expression.H"
25 : #include "Numeric/Stencil.H"
26 :
27 : namespace Integrator
28 : {
29 : class Dendrite : virtual public Integrator
30 : {
31 : public:
32 : static constexpr const char* name = "dendrite";
33 :
34 : // Empty constructor
35 1 : Dendrite() : Integrator()
36 1 : {}
37 :
38 : // Constructor that triggers parse
39 1 : Dendrite(IO::ParmParse& pp) : Dendrite()
40 : {
41 1 : Parse(*this, pp);
42 1 : }
43 :
44 2 : ~Dendrite()
45 1 : {
46 1 : delete ic_temp;
47 1 : delete ic_phi;
48 1 : delete bc_temp;
49 1 : delete bc_phi;
50 2 : }
51 :
52 1 : static void Parse(Dendrite& value, IO::ParmParse& pp)
53 : {
54 6 : pp_query_required("alpha", value.alpha); // Pre-multiplier of "m" barrier height
55 6 : pp_query_required("delta", value.delta); // Anisotropy factor
56 6 : pp_query_required("gamma", value.gamma); // Anisotropic temperature coupling factor
57 6 : pp_query_default("diffusion", value.diffusion,1.0); // Thermal constant
58 :
59 6 : pp_query_required("eps", value.eps); // Diffuse boundary width
60 6 : pp_query_required("tau", value.tau); // Diffusive timescale
61 :
62 5 : pp_query_default("theta", value.theta, 0.0); // Orientation about z axis (Deg)
63 :
64 1 : Eigen::Matrix3d R3d;
65 1 : R3d = Eigen::AngleAxisd(value.theta*Set::Constant::Pi/180.0, Eigen::Vector3d::UnitZ());
66 1 : value.R = Set::reduce(R3d);
67 :
68 : // Refinement criteria for temperature
69 6 : pp_query_default("heat.refinement_threshold", value.refinement_threshold_temp, 0.01);
70 : // Refinement criteria for phi
71 5 : pp_query_default("phi.refinement_threshold", value.refinement_threshold_phi, 0.01);
72 :
73 : // Expressions for ICs
74 2 : value.ic_temp = new IC::Expression(value.geom, pp, "ic.temp");
75 2 : value.ic_phi = new IC::Expression(value.geom, pp, "ic.phi");
76 :
77 : // boundary conditions for temperature field
78 2 : pp.select_default<BC::Constant>("bc.temp", value.bc_temp, 1);
79 : // boundary conditions for :math:`\phi` field
80 2 : pp.select_default<BC::Constant>("bc.phi", value.bc_phi, 1);
81 :
82 3 : value.RegisterNewFab(value.temp_mf, value.bc_temp, value.number_of_components, value.number_of_ghost_cells, "Temp", true);
83 3 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, value.number_of_components, value.number_of_ghost_cells, "Temp_old", false);
84 3 : value.RegisterNewFab(value.phi_mf, value.bc_phi, value.number_of_components, value.number_of_ghost_cells, "phi", true);
85 3 : value.RegisterNewFab(value.phi_old_mf, value.bc_phi, value.number_of_components, value.number_of_ghost_cells, "phi_old", false);
86 1 : }
87 : protected:
88 :
89 : // Use the ic object to initialize the temperature field
90 4 : void Initialize(int lev)
91 : {
92 4 : ic_temp->Initialize(lev, temp_old_mf);
93 4 : ic_phi->Initialize(lev, phi_old_mf);
94 4 : ic_temp->Initialize(lev, temp_mf);
95 4 : ic_phi->Initialize(lev, phi_mf);
96 4 : }
97 :
98 :
99 : // Integrate the heat equation
100 14 : void Advance(int lev, Set::Scalar /*time*/, Set::Scalar dt)
101 : {
102 : // Swap the old temp fab and the new temp fab so we use
103 : // the new one.
104 14 : std::swap(*temp_mf[lev], *temp_old_mf[lev]);
105 14 : std::swap(*phi_mf[lev], *phi_old_mf[lev]);
106 :
107 : // Get the cell size corresponding to this level
108 14 : const Set::Scalar* DX = geom[lev].CellSize();
109 :
110 : // Iterate over all of the patches on this level
111 88 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
112 : {
113 : // Get the box (index dimensions) for this patch
114 74 : const amrex::Box& bx = mfi.tilebox();
115 :
116 : // Get an array-accessible handle to the data on this patch.
117 74 : amrex::Array4<const Set::Scalar> const& temp = (*temp_old_mf[lev]).array(mfi);
118 74 : amrex::Array4<Set::Scalar> const& temp_new = (*temp_mf[lev]).array(mfi);
119 74 : amrex::Array4<const Set::Scalar> const& phi = (*phi_old_mf[lev]).array(mfi);
120 74 : amrex::Array4<Set::Scalar> const& phi_new = (*phi_mf[lev]).array(mfi);
121 :
122 : // Iterate over the grid on this patch
123 74 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
124 : {
125 :
126 10752 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
127 :
128 10752 : gradphi = R * gradphi;
129 :
130 : //
131 : // Calculate anisotropy function
132 : // [ equation 2.12 in reference]
133 : //
134 10752 : Set::Scalar v44 = AMREX_D_TERM(
135 : gradphi(0) * gradphi(0) * gradphi(0) * gradphi(0),
136 : +gradphi(1) * gradphi(1) * gradphi(1) * gradphi(1),
137 : +gradphi(2) * gradphi(2) * gradphi(2) * gradphi(2));
138 10752 : Set::Scalar v22 = gradphi.squaredNorm(); v22 *= v22;
139 10752 : Set::Scalar sigma = 1.0 - delta * (1.0 - v44 / (v22 + 1E-12));
140 :
141 : //
142 : // Calculate anisotropic barrier height
143 : // [ unnumbered equation on page 67 in the reference ]
144 : //
145 10752 : Set::Scalar m = -(alpha / Set::Constant::Pi) * std::atan(gamma * sigma * temp(i, j, k));
146 :
147 :
148 : //
149 : // Evolve order paramter
150 : // [ equation 2.10 in the reference ]
151 : //
152 10752 : Set::Scalar lapphi = Numeric::Laplacian(phi, i, j, k, 0, DX);
153 32256 : Set::Scalar Dphi = (eps * eps * lapphi + phi(i, j, k) * (1.0 - phi(i, j, k)) * (phi(i, j, k) - 0.5 + m)) / tau;
154 21504 : phi_new(i, j, k) = phi(i, j, k) + dt * Dphi;
155 :
156 : //
157 : // Evolve temperature
158 : // [ equation 2.11 in the reference ]
159 : //
160 10752 : Set::Scalar laptemp = Numeric::Laplacian(temp, i, j, k, 0, DX);
161 21504 : temp_new(i, j, k) = temp(i, j, k) + dt * (diffusion*laptemp + Dphi);
162 10752 : });
163 14 : }
164 14 : }
165 :
166 : // Tag cells for mesh refinement based on temperature gradient
167 9 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
168 : {
169 : // Get cell dimensions as done above.
170 9 : const Set::Scalar* DX = geom[lev].CellSize();
171 : // Calculate the diagonal.
172 9 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
173 :
174 : // Iterate over the patches on this level
175 33 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
176 : {
177 : // Get the box and handles as done above.
178 24 : const amrex::Box& bx = mfi.tilebox();
179 24 : amrex::Array4<char> const& tags = a_tags.array(mfi);
180 24 : amrex::Array4<Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
181 24 : amrex::Array4<Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
182 :
183 : // Iterate over the grid as done above.
184 24 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
185 : {
186 : // Calculate the temperature gradient.
187 17024 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
188 17024 : Set::Vector grad_phi = Numeric::Gradient(phi, i, j, k, 0, DX);
189 :
190 : // Is the gradient * cell_size too big? If so, then
191 : // mark this cell as needing refinement.
192 17024 : if (grad_temp.lpNorm<2>() * dr > refinement_threshold_temp)
193 278 : tags(i, j, k) = amrex::TagBox::SET;
194 17024 : if (grad_phi.lpNorm<2>() * dr > refinement_threshold_phi)
195 278 : tags(i, j, k) = amrex::TagBox::SET;
196 17024 : });
197 9 : }
198 9 : }
199 :
200 : protected:
201 : Set::Field<Set::Scalar> temp_mf;
202 : Set::Field<Set::Scalar> temp_old_mf;
203 : Set::Field<Set::Scalar> phi_mf;
204 : Set::Field<Set::Scalar> phi_old_mf;
205 :
206 : Set::Field<Set::Scalar> &eta_mf = phi_mf;
207 : Set::Field<Set::Scalar> &eta_old_mf = phi_old_mf;
208 :
209 : private:
210 : int number_of_components = 1; // Number of components
211 : int number_of_ghost_cells = 1; // Number of ghost cells
212 :
213 : // calculation of m
214 : Set::Scalar alpha, delta, gamma, diffusion=NAN;
215 : // evolution of p
216 : Set::Scalar eps, tau;
217 : // orientation
218 : Set::Scalar theta = NAN;
219 : Set::Matrix R = Set::Matrix::Identity();
220 :
221 : Set::Scalar refinement_threshold_temp = NAN; // Criterion for cell refinement
222 : Set::Scalar refinement_threshold_phi = NAN; // Criterion for cell refinement
223 :
224 : IC::IC<Set::Scalar>* ic_temp, * ic_phi;
225 : BC::BC<Set::Scalar>* bc_temp, * bc_phi;
226 : };
227 : } // namespace Integrator
228 : #endif
|
