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