Line data Source code
1 : //
2 : // This is a simple implementation of an Allen-Cahn equation governed by
3 : //
4 : // .. math::
5 : //
6 : // \frac{\partial\alpha}{\partial t} =
7 : // - L \Big(\frac{\lambda}{\epsilon}(2.0\alpha
8 : // - 6\alpha^2 + 4\alpha^3) + \epsilon\,\kappa\,\Delta \alpha \Big)
9 : //
10 : // where :math:`\alpha(x,t)` is the order parameter.
11 : //
12 :
13 : #ifndef INTEGRATOR_ALLENCAHN_H // Include guards
14 : #define INTEGRATOR_ALLENCAHN_H //
15 :
16 : // Standard library includes
17 : #include <iostream>
18 : #include <fstream>
19 : #include <iomanip>
20 :
21 : // AMReX Includes
22 : #include "AMReX.H"
23 : #include "AMReX_ParallelDescriptor.H"
24 : #include "AMReX_ParmParse.H"
25 :
26 : // Alamo Includes
27 : #include "IO/ParmParse.H"
28 : #include "Integrator/Integrator.H"
29 : #include "BC/Constant.H"
30 : #include "IC/IC.H"
31 : #include "IC/Sphere.H"
32 : #include "IC/Constant.H"
33 : #include "IC/Expression.H"
34 : #include "IC/BMP.H"
35 : #include "IC/PNG.H"
36 : #include "IC/PSRead.H"
37 : #include "Numeric/Stencil.H"
38 : #include "IC/Random.H"
39 : #include "Unit/Unit.H"
40 :
41 : namespace Integrator
42 : {
43 :
44 : /// This is the definition of the Allen Cahn class
45 : class AllenCahn : virtual public Integrator
46 : {
47 : public:
48 : static constexpr const char* name = "allencahn";
49 :
50 : // Empty constructor
51 3 : AllenCahn() : Integrator()
52 3 : {}
53 :
54 : // Constructor that triggers parse
55 3 : AllenCahn(IO::ParmParse& pp) : AllenCahn()
56 : {
57 3 : Parse(*this, pp);
58 3 : }
59 :
60 6 : ~AllenCahn()
61 3 : {
62 3 : delete ic;
63 3 : delete bc;
64 6 : }
65 :
66 : // The Parse function initializes the AllenCahn object using
67 : // a parser, pp.
68 : // Note that this is a static function, which means it does not have
69 : // direct access to member variables. Instead, it initializes the variables
70 : // inside the argument, "value", and so all references to member items are
71 : // prefixed by "value."
72 3 : static void Parse(AllenCahn& value, IO::ParmParse& pp)
73 : {
74 : // Criterion for mesh refinement [0.01]
75 6 : pp_query_default("refinement_threshold", value.refinement_threshold, 0.01);
76 :
77 : // Value for :math:`L` (mobility)
78 12 : pp.query_default("ch.L",value.ch.L, "1.0", 1/ Unit::Pressure()/Unit::Time());
79 : // Value for :math:`\epsilon` (diffuse boundary width)
80 : // pp.query_default("ch.eps",value.ch.eps, "0.1", Unit::Length());
81 : // Value for :math:`\kappa` (Interface energy parameter)
82 12 : pp.query_default("ch.grad",value.ch.kappa, "1.0", Unit::Energy()/Unit::Length());
83 : // Value for :math:`\lambda` (Chemical potential coefficient)
84 12 : pp_query_default("ch.chempot",value.ch.lambda, "1.0", Unit::Pressure());
85 : // Force directional growth: 0=no growth, 1=only positive, -1=only negative
86 12 : pp_query_validate("ch.direction",value.ch.direction, {0,1,-1});
87 : // Time to start forcing directional growth
88 9 : pp_query_default("ch.direction_tstart",value.ch.direction_tstart, 0.0);
89 :
90 :
91 15 : std::pair<std::string,Set::Scalar> ch_coeffs[2];
92 12 : pp.query_exactly<2>({"ch.lambda","ch.kappa","ch.sigma","ch.eps"}, ch_coeffs,
93 3 : {Unit::Pressure(), Unit::Energy()/Unit::Length(), Unit::Energy()/Unit::Area(), Unit::Length()});
94 :
95 3 : if (ch_coeffs[0].first == "ch.lambda" && ch_coeffs[1].first == "ch.kappa")
96 : {
97 2 : value.ch.lambda = ch_coeffs[0].second;
98 2 : value.ch.kappa = ch_coeffs[1].second;
99 : }
100 1 : else if (ch_coeffs[0].first == "ch.sigma" && ch_coeffs[1].first == "ch.eps")
101 : {
102 1 : Set::Scalar sigma = ch_coeffs[0].second, eps = ch_coeffs[1].second;
103 1 : value.ch.lambda = 4.0 * sigma / eps;
104 1 : value.ch.kappa = sigma * eps / 2.0;
105 : }
106 :
107 :
108 : // Set the initial condition for the alpha field
109 6 : pp.select_default<IC::Sphere, IC::Constant, IC::Expression,IC::BMP, IC::PNG, IC::Random, IC::PSRead>("alpha.ic", value.ic, value.geom);
110 :
111 :
112 : // Use a constant BC object for temperature
113 : // value.bc = new BC::Constant(1);
114 : // :ref:`BC::Constant` parameters
115 6 : pp.select_default<BC::Constant>("alpha.bc", value.bc, 1);
116 :
117 : // Register the temperature and old temperature fields.
118 : // alpha_mf and alpha_old_mf are defined near the bottom of this Header file.
119 9 : value.RegisterNewFab(value.alpha_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "alpha", true);
120 9 : value.RegisterNewFab(value.alpha_old_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "alpha_old", false);
121 12 : }
122 :
123 : protected:
124 :
125 : // Use the ic object to initialize the order parameter
126 10 : void Initialize(int lev)
127 : {
128 10 : ic->Initialize(lev, alpha_mf);
129 10 : ic->Initialize(lev, alpha_old_mf);
130 10 : }
131 :
132 : // Integrate the Allen Cahn equation
133 82700 : void Advance(int lev, Set::Scalar time, Set::Scalar dt)
134 : {
135 82700 : std::swap(*alpha_mf[lev], *alpha_old_mf[lev]);
136 :
137 82700 : const Set::Scalar *DX = this->geom[lev].CellSize();
138 : // Evolve alpha
139 165400 : for (amrex::MFIter mfi(*alpha_mf[lev], true); mfi.isValid(); ++mfi)
140 : {
141 82700 : const amrex::Box& bx = mfi.tilebox();
142 82700 : Set::Patch<Set::Scalar> alpha_new = alpha_mf.Patch(lev,mfi);
143 82700 : Set::Patch<const Set::Scalar> alpha = alpha_old_mf.Patch(lev,mfi);
144 :
145 82700 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
146 : {
147 265411456 : Set::Scalar driving_force = 0.0;
148 :
149 : // Chemical potential
150 530822912 : Set::Scalar alpha_2 = alpha(i, j, k) * alpha(i, j, k);
151 265411456 : Set::Scalar alpha_3 = alpha_2 * alpha(i, j, k);
152 265411456 : driving_force += ch.lambda * (2.0 * alpha(i, j, k) - 6.0 * alpha_2 + 4.0 * alpha_3);
153 :
154 : // Gradient
155 265411456 : Set::Scalar alpha_lap = Numeric::Laplacian(alpha, i, j, k, 0, DX);
156 265411456 : driving_force -= ch.kappa * alpha_lap;
157 :
158 265411456 : if (time >= ch.direction_tstart)
159 : {
160 265411456 : if (ch.direction == 1)
161 0 : driving_force = std::min(0.0, driving_force);
162 265411456 : else if (ch.direction == -1)
163 0 : driving_force = std::max(0.0, driving_force);
164 : }
165 :
166 : // Update
167 530822912 : alpha_new(i, j, k) =
168 530822912 : alpha(i, j, k) - ch.L * dt * driving_force;
169 265411456 : });
170 82700 : }
171 82700 : }
172 :
173 : // Tag cells for mesh refinement based on temperature gradient
174 29714 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
175 : {
176 : // Get cell dimensions as done above.
177 29714 : const Set::Scalar* DX = geom[lev].CellSize();
178 : // Calculate the diagonal.
179 29714 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
180 :
181 : // Iterate over the patches on this level
182 59428 : for (amrex::MFIter mfi(*alpha_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
183 : {
184 : // Get the box and handles as done above.
185 29714 : const amrex::Box& bx = mfi.tilebox();
186 29714 : amrex::Array4<char> const& tags = a_tags.array(mfi);
187 29714 : amrex::Array4<Set::Scalar> const& alpha = (*alpha_mf[lev]).array(mfi);
188 :
189 : // Iterate over the grid as done above.
190 29714 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
191 : {
192 : // Calculate the temperature gradient.
193 28848736 : Set::Vector grad = Numeric::Gradient(alpha, i, j, k, 0, DX);
194 :
195 : // Is the gradient * cell_size too big? If so, then
196 : // mark this cell as needing refinement.
197 28848736 : if (grad.lpNorm<2>() * dr > refinement_threshold)
198 41137100 : tags(i, j, k) = amrex::TagBox::SET;
199 28848736 : });
200 29714 : }
201 29714 : }
202 :
203 : protected:
204 : Set::Field<Set::Scalar> alpha_mf; /// Temperature field variable (current timestep)
205 : Set::Field<Set::Scalar> alpha_old_mf; /// Temperature field variable (previous timestep)
206 : IC::IC<Set::Scalar>* ic; /// Object used to initialize temperature field
207 :
208 : Set::Field<Set::Scalar> &eta_mf = alpha_mf;
209 : Set::Field<Set::Scalar> &eta_old_mf = alpha_old_mf;
210 :
211 : struct {
212 : Set::Scalar L = NAN;
213 : Set::Scalar eps = NAN;
214 : Set::Scalar kappa = NAN;
215 : Set::Scalar lambda = NAN;
216 : Set::Scalar direction_tstart = NAN;
217 : int direction = 0;
218 : } ch;
219 :
220 : private:
221 : int number_of_components = 1; // Number of components
222 : int number_of_ghost_cells = 2; // Number of ghost cells
223 :
224 : Set::Scalar refinement_threshold = NAN; // Criterion for cell refinement
225 :
226 : BC::BC<Set::Scalar>* bc; // Object used to update temp field boundary ghost cells
227 : };
228 : } // namespace Integrator
229 : #endif
|
