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 :
40 : namespace Integrator
41 : {
42 :
43 : /// This is the definition of the Allen Cahn class
44 : class AllenCahn : virtual public Integrator
45 : {
46 : public:
47 : static constexpr const char* name = "allencahn";
48 :
49 : // Empty constructor
50 1 : AllenCahn() : Integrator()
51 1 : {}
52 :
53 : // Constructor that triggers parse
54 1 : AllenCahn(IO::ParmParse& pp) : AllenCahn()
55 : {
56 1 : Parse(*this, pp);
57 1 : }
58 :
59 2 : ~AllenCahn()
60 1 : {
61 1 : delete ic;
62 1 : delete bc;
63 2 : }
64 :
65 : // The Parse function initializes the AllenCahn object using
66 : // a parser, pp.
67 : // Note that this is a static function, which means it does not have
68 : // direct access to member variables. Instead, it initializes the variables
69 : // inside the argument, "value", and so all references to member items are
70 : // prefixed by "value."
71 1 : static void Parse(AllenCahn& value, IO::ParmParse& pp)
72 : {
73 : // Criterion for mesh refinement [0.01]
74 6 : pp_query_default("refinement_threshold", value.refinement_threshold, 0.01);
75 :
76 : // Value for :math:`L` (mobility)
77 6 : pp_query_default("ch.L",value.ch.L, 1.0);
78 : // Value for :math:`\epsilon` (diffuse boundary width)
79 6 : pp_query_default("ch.eps",value.ch.eps, 0.1);
80 : // Value for :math:`\kappa` (Interface energy parameter)
81 6 : pp_query_default("ch.grad",value.ch.grad, 1.0);
82 : // Value for :math:`\lambda` (Chemical potential coefficient)
83 6 : pp_query_default("ch.chempot",value.ch.chempot, 1.0);
84 : // Force directional growth: 0=no growth, 1=only positive, -1=only negative
85 8 : pp_query_validate("ch.direction",value.ch.direction, {0,1,-1});
86 : // Time to start forcing directional growth
87 5 : pp_query_default("ch.direction_tstart",value.ch.direction_tstart, 0.0);
88 :
89 : // Set the initial condition for the alpha field
90 2 : pp.select_default<IC::Sphere, IC::Constant, IC::Expression, IC::BMP, IC::PNG, IC::Random, IC::PSRead>("alpha.ic", value.ic, value.geom);
91 :
92 :
93 : // Use a constant BC object for temperature
94 : // value.bc = new BC::Constant(1);
95 : // :ref:`BC::Constant` parameters
96 2 : pp.select_default<BC::Constant>("alpha.bc", value.bc, 1);
97 :
98 : // Register the temperature and old temperature fields.
99 : // alpha_mf and alpha_old_mf are defined near the bottom of this Header file.
100 3 : value.RegisterNewFab(value.alpha_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "alpha", true);
101 3 : value.RegisterNewFab(value.alpha_old_mf, value.bc, value.number_of_components, value.number_of_ghost_cells, "alpha_old", false);
102 1 : }
103 :
104 : protected:
105 :
106 : // Use the ic object to initialize the order parameter
107 3 : void Initialize(int lev)
108 : {
109 3 : ic->Initialize(lev, alpha_mf);
110 3 : ic->Initialize(lev, alpha_old_mf);
111 3 : }
112 :
113 : // Integrate the Allen Cahn equation
114 700 : void Advance(int lev, Set::Scalar time, Set::Scalar dt)
115 : {
116 700 : std::swap(*alpha_mf[lev], *alpha_old_mf[lev]);
117 :
118 700 : const Set::Scalar *DX = this->geom[lev].CellSize();
119 : // Evolve alpha
120 1400 : for (amrex::MFIter mfi(*alpha_mf[lev], true); mfi.isValid(); ++mfi)
121 : {
122 700 : const amrex::Box& bx = mfi.tilebox();
123 700 : Set::Patch<Set::Scalar> alpha_new = alpha_mf.Patch(lev,mfi);
124 700 : Set::Patch<const Set::Scalar> alpha = alpha_old_mf.Patch(lev,mfi);
125 :
126 700 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
127 : {
128 19976192 : Set::Scalar driving_force = 0.0;
129 : // Chemical potential
130 39952384 : Set::Scalar alpha_2 = alpha(i, j, k) * alpha(i, j, k);
131 19976192 : Set::Scalar alpha_3 = alpha_2 * alpha(i, j, k);
132 19976192 : driving_force += (ch.chempot / ch.eps) * (2.0 * alpha(i, j, k) - 6.0 * alpha_2 + 4.0 * alpha_3);
133 : // Gradient
134 19976192 : Set::Scalar alpha_lap = Numeric::Laplacian(alpha, i, j, k, 0, DX);
135 19976192 : driving_force -= ch.eps * ch.grad * alpha_lap;
136 :
137 19976192 : if (time >= ch.direction_tstart)
138 : {
139 19976192 : if (ch.direction == 1)
140 0 : driving_force = std::min(0.0, driving_force);
141 19976192 : else if (ch.direction == -1)
142 0 : driving_force = std::max(0.0, driving_force);
143 : }
144 :
145 : // Update
146 39952384 : alpha_new(i, j, k) =
147 39952384 : alpha(i, j, k) - ch.L * dt * driving_force;
148 19976192 : });
149 700 : }
150 700 : }
151 :
152 : // Tag cells for mesh refinement based on temperature gradient
153 204 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
154 : {
155 : // Get cell dimensions as done above.
156 204 : const Set::Scalar* DX = geom[lev].CellSize();
157 : // Calculate the diagonal.
158 204 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
159 :
160 : // Iterate over the patches on this level
161 408 : for (amrex::MFIter mfi(*alpha_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
162 : {
163 : // Get the box and handles as done above.
164 204 : const amrex::Box& bx = mfi.tilebox();
165 204 : amrex::Array4<char> const& tags = a_tags.array(mfi);
166 204 : amrex::Array4<Set::Scalar> const& alpha = (*alpha_mf[lev]).array(mfi);
167 :
168 : // Iterate over the grid as done above.
169 204 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
170 : {
171 : // Calculate the temperature gradient.
172 2028544 : Set::Vector grad = Numeric::Gradient(alpha, i, j, k, 0, DX);
173 :
174 : // Is the gradient * cell_size too big? If so, then
175 : // mark this cell as needing refinement.
176 2028544 : if (grad.lpNorm<2>() * dr > refinement_threshold)
177 3425626 : tags(i, j, k) = amrex::TagBox::SET;
178 2028544 : });
179 204 : }
180 204 : }
181 :
182 : protected:
183 : Set::Field<Set::Scalar> alpha_mf; /// Temperature field variable (current timestep)
184 : Set::Field<Set::Scalar> alpha_old_mf; /// Temperature field variable (previous timestep)
185 : IC::IC<Set::Scalar>* ic; /// Object used to initialize temperature field
186 :
187 : Set::Field<Set::Scalar> &eta_mf = alpha_mf;
188 : Set::Field<Set::Scalar> &eta_old_mf = alpha_old_mf;
189 :
190 : struct {
191 : Set::Scalar L = NAN;
192 : Set::Scalar eps = NAN;
193 : Set::Scalar grad = NAN;
194 : Set::Scalar chempot = NAN;
195 : Set::Scalar direction_tstart = NAN;
196 : int direction = 0;
197 : } ch;
198 :
199 : private:
200 : int number_of_components = 1; // Number of components
201 : int number_of_ghost_cells = 2; // Number of ghost cells
202 :
203 : Set::Scalar refinement_threshold = NAN; // Criterion for cell refinement
204 :
205 : BC::BC<Set::Scalar>* bc; // Object used to update temp field boundary ghost cells
206 : };
207 : } // namespace Integrator
208 : #endif
|
