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