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