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 : // Empty constructor
33 2 : Dendrite() : Integrator()
34 2 : {}
35 :
36 : // Constructor that triggers parse
37 2 : Dendrite(IO::ParmParse& pp) : Dendrite()
38 : {
39 2 : Parse(*this, pp);
40 2 : }
41 :
42 2 : ~Dendrite()
43 1 : {
44 1 : delete ic_temp;
45 1 : delete ic_phi;
46 1 : delete bc_temp;
47 1 : delete bc_phi;
48 2 : }
49 :
50 2 : static void Parse(Dendrite& value, IO::ParmParse& pp)
51 : {
52 12 : pp_query_required("alpha", value.alpha); // Pre-multiplier of "m" barrier height
53 12 : pp_query_required("delta", value.delta); // Anisotropy factor
54 12 : pp_query_required("gamma", value.gamma); // Anisotropic temperature coupling factor
55 :
56 12 : pp_query_required("eps", value.eps); // Diffuse boundary width
57 12 : pp_query_required("tau", value.tau); // Diffusive timescale
58 :
59 : // Refinement criteria for temperature
60 12 : pp_query_default("heat.refinement_threshold", value.refinement_threshold_temp, 0.01);
61 : // Refinement criteria for phi
62 10 : pp_query_default("phi.refinement_threshold", value.refinement_threshold_phi, 0.01);
63 :
64 : // Expressions for ICs
65 4 : value.ic_temp = new IC::Expression(value.geom, pp, "ic.temp");
66 4 : value.ic_phi = new IC::Expression(value.geom, pp, "ic.phi");
67 :
68 : // Use a constant BC object for temperature
69 2 : value.bc_temp = new BC::Constant(1);
70 2 : value.bc_phi = new BC::Constant(1);
71 12 : pp_queryclass("bc.temp", *static_cast<BC::Constant*>(value.bc_temp));
72 10 : pp_queryclass("bc.phi", *static_cast<BC::Constant*>(value.bc_phi));
73 :
74 6 : value.RegisterNewFab(value.temp_mf, value.bc_temp, value.number_of_components, value.number_of_ghost_cells, "Temp", true);
75 6 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, value.number_of_components, value.number_of_ghost_cells, "Temp_old", false);
76 6 : value.RegisterNewFab(value.phi_mf, value.bc_phi, value.number_of_components, value.number_of_ghost_cells, "phi", true);
77 6 : value.RegisterNewFab(value.phi_old_mf, value.bc_phi, value.number_of_components, value.number_of_ghost_cells, "phi_old", false);
78 2 : }
79 : protected:
80 :
81 : // Use the ic object to initialize the temperature field
82 8 : void Initialize(int lev)
83 : {
84 8 : ic_temp->Initialize(lev, temp_old_mf);
85 8 : ic_phi->Initialize(lev, phi_old_mf);
86 8 : ic_temp->Initialize(lev, temp_mf);
87 8 : ic_phi->Initialize(lev, phi_mf);
88 8 : }
89 :
90 :
91 : // Integrate the heat equation
92 18332 : void Advance(int lev, Set::Scalar /*time*/, Set::Scalar dt)
93 : {
94 : // Swap the old temp fab and the new temp fab so we use
95 : // the new one.
96 18332 : std::swap(*temp_mf[lev], *temp_old_mf[lev]);
97 18332 : std::swap(*phi_mf[lev], *phi_old_mf[lev]);
98 :
99 : // Get the cell size corresponding to this level
100 18332 : const Set::Scalar* DX = geom[lev].CellSize();
101 :
102 : // Iterate over all of the patches on this level
103 248783 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
104 : {
105 : // Get the box (index dimensions) for this patch
106 230452 : const amrex::Box& bx = mfi.tilebox();
107 :
108 : // Get an array-accessible handle to the data on this patch.
109 230452 : amrex::Array4<const Set::Scalar> const& temp = (*temp_old_mf[lev]).array(mfi);
110 230452 : amrex::Array4<Set::Scalar> const& temp_new = (*temp_mf[lev]).array(mfi);
111 230452 : amrex::Array4<const Set::Scalar> const& phi = (*phi_old_mf[lev]).array(mfi);
112 230452 : amrex::Array4<Set::Scalar> const& phi_new = (*phi_mf[lev]).array(mfi);
113 :
114 : // Iterate over the grid on this patch
115 230452 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
116 : {
117 :
118 39560838 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
119 :
120 : //
121 : // Calculate anisotropy function
122 : // [ equation 2.12 in reference]
123 : //
124 39560837 : Set::Scalar v44 = AMREX_D_TERM(
125 : gradphi(0) * gradphi(0) * gradphi(0) * gradphi(0),
126 : +gradphi(1) * gradphi(1) * gradphi(1) * gradphi(1),
127 : +gradphi(2) * gradphi(2) * gradphi(2) * gradphi(2));
128 39560837 : Set::Scalar v22 = gradphi.squaredNorm(); v22 *= v22;
129 39560837 : Set::Scalar sigma = 1.0 - delta * (1.0 - v44 / (v22 + 1E-12));
130 :
131 : //
132 : // Calculate anisotropic barrier height
133 : // [ unnumbered equation on page 67 in the reference ]
134 : //
135 79121674 : Set::Scalar m = -(alpha / Set::Constant::Pi) * std::atan(gamma * sigma * temp(i, j, k));
136 :
137 :
138 : //
139 : // Evolve order paramter
140 : // [ equation 2.10 in the reference ]
141 : //
142 39560837 : Set::Scalar lapphi = Numeric::Laplacian(phi, i, j, k, 0, DX);
143 118682511 : Set::Scalar Dphi = (eps * eps * lapphi + phi(i, j, k) * (1.0 - phi(i, j, k)) * (phi(i, j, k) - 0.5 + m)) / tau;
144 79121674 : phi_new(i, j, k) = phi(i, j, k) + dt * Dphi;
145 :
146 : //
147 : // Evolve temperature
148 : // [ equation 2.11 in the reference ]
149 : //
150 39560837 : Set::Scalar laptemp = Numeric::Laplacian(temp, i, j, k, 0, DX);
151 79121674 : temp_new(i, j, k) = temp(i, j, k) + dt * (laptemp + Dphi);
152 39560837 : });
153 18331 : }
154 18331 : }
155 :
156 : // Tag cells for mesh refinement based on temperature gradient
157 1063 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
158 : {
159 : // Get cell dimensions as done above.
160 1063 : const Set::Scalar* DX = geom[lev].CellSize();
161 : // Calculate the diagonal.
162 1063 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
163 :
164 : // Iterate over the patches on this level
165 8090 : for (amrex::MFIter mfi(*temp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
166 : {
167 : // Get the box and handles as done above.
168 7027 : const amrex::Box& bx = mfi.tilebox();
169 7027 : amrex::Array4<char> const& tags = a_tags.array(mfi);
170 7027 : amrex::Array4<Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
171 7027 : amrex::Array4<Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
172 :
173 : // Iterate over the grid as done above.
174 7027 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
175 : {
176 : // Calculate the temperature gradient.
177 1848704 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
178 1848704 : Set::Vector grad_phi = Numeric::Gradient(phi, i, j, k, 0, DX);
179 :
180 : // Is the gradient * cell_size too big? If so, then
181 : // mark this cell as needing refinement.
182 1848704 : if (grad_temp.lpNorm<2>() * dr > refinement_threshold_temp)
183 716616 : tags(i, j, k) = amrex::TagBox::SET;
184 1848704 : if (grad_phi.lpNorm<2>() * dr > refinement_threshold_phi)
185 246940 : tags(i, j, k) = amrex::TagBox::SET;
186 1848704 : });
187 1063 : }
188 1063 : }
189 :
190 : protected:
191 : Set::Field<Set::Scalar> temp_mf;
192 : Set::Field<Set::Scalar> temp_old_mf;
193 : Set::Field<Set::Scalar> phi_mf;
194 : Set::Field<Set::Scalar> phi_old_mf;
195 :
196 : private:
197 : int number_of_components = 1; // Number of components
198 : int number_of_ghost_cells = 1; // Number of ghost cells
199 :
200 :
201 : // calculation of m
202 : Set::Scalar alpha, delta, gamma;
203 : // evolution of p
204 : Set::Scalar eps, tau;
205 :
206 : Set::Scalar refinement_threshold_temp = NAN; // Criterion for cell refinement
207 : Set::Scalar refinement_threshold_phi = NAN; // Criterion for cell refinement
208 :
209 : IC::IC<Set::Scalar>* ic_temp, * ic_phi;
210 : BC::BC<Set::Scalar>* bc_temp, * bc_phi;
211 : };
212 : } // namespace Integrator
213 : #endif
|
