Line data Source code
1 : #include <AMReX_MLPoisson.H>
2 :
3 : #ifdef ALAMO_FFT
4 : #include <AMReX_FFT.H>
5 : #endif
6 :
7 : #include "CahnHilliard.H"
8 : #include "BC/Constant.H"
9 : #include "IO/ParmParse.H"
10 : #include "IC/Random.H"
11 : #include "Numeric/Stencil.H"
12 : #include "Set/Set.H"
13 :
14 : namespace Integrator
15 : {
16 0 : CahnHilliard::CahnHilliard() : Integrator()
17 : {
18 0 : }
19 0 : CahnHilliard::~CahnHilliard()
20 : {
21 0 : delete ic;
22 0 : delete bc;
23 0 : }
24 :
25 0 : void CahnHilliard::Parse(CahnHilliard &value, IO::ParmParse &pp)
26 : {
27 : // Interface energy
28 0 : pp.query_default("gamma",value.gamma, 0.0005);
29 : // Mobility
30 0 : pp.query_default("L", value.L, 1.0);
31 : // Regridding criterion
32 0 : pp.query_default("refinement_threshold",value.refinement_threshold, 1E100);
33 :
34 : // initial condition for :math:`\eta`
35 0 : pp.select_default<IC::Random>("eta.ic", value.ic, value.geom);
36 : // boundary condition for :math:`\eta`
37 0 : pp.select_default<BC::Constant>("eta.bc", value.bc, 1);
38 :
39 : // Which method to use - realspace or spectral method.
40 0 : pp.query_validate("method",value.method,{"realspace","spectral"});
41 :
42 0 : value.RegisterNewFab(value.etanew_mf, value.bc, 1, 1, "eta",true);
43 0 : value.RegisterNewFab(value.intermediate, value.bc, 1, 1, "int",true);
44 :
45 0 : if (value.method == "realspace")
46 0 : value.RegisterNewFab(value.etaold_mf, value.bc, 1, 1, "eta_old",false);
47 0 : }
48 :
49 : void
50 0 : CahnHilliard::Advance(int lev, Set::Scalar time, Set::Scalar dt)
51 : {
52 0 : if (method == "realspace")
53 0 : AdvanceReal(lev, time, dt);
54 0 : else if (method == "spectral")
55 0 : AdvanceSpectral(lev, time, dt);
56 : else
57 0 : Util::Abort(INFO,"Invalid method: ",method);
58 0 : }
59 :
60 :
61 : void
62 0 : CahnHilliard::AdvanceReal (int lev, Set::Scalar /*time*/, Set::Scalar dt)
63 : {
64 0 : std::swap(etaold_mf[lev], etanew_mf[lev]);
65 0 : const Set::Scalar* DX = geom[lev].CellSize();
66 0 : for ( amrex::MFIter mfi(*etanew_mf[lev],true); mfi.isValid(); ++mfi )
67 : {
68 0 : const amrex::Box& bx = mfi.tilebox();
69 0 : amrex::Array4<const amrex::Real> const& eta = etaold_mf[lev]->array(mfi);
70 0 : amrex::Array4<amrex::Real> const& inter = intermediate[lev]->array(mfi);
71 0 : amrex::Array4<amrex::Real> const& etanew = etanew_mf[lev]->array(mfi);
72 :
73 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
74 : {
75 0 : Set::Scalar lap_eta = Numeric::Laplacian(eta,i,j,k,0,DX);
76 :
77 :
78 0 : inter(i,j,k) =
79 0 : eta(i,j,k)*eta(i,j,k)*eta(i,j,k)
80 0 : - eta(i,j,k)
81 0 : - gamma*lap_eta;
82 :
83 :
84 0 : etanew(i,j,k) = eta(i,j,k) - dt*inter(i,j,k); // Allen Cahn
85 0 : });
86 :
87 0 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k){
88 0 : Set::Scalar lap_inter = Numeric::Laplacian(inter,i,j,k,0,DX);
89 :
90 0 : etanew(i,j,k) = eta(i,j,k) + dt*lap_inter;
91 0 : });
92 0 : }
93 0 : }
94 :
95 : #ifdef ALAMO_FFT
96 : void
97 : CahnHilliard::AdvanceSpectral (int lev, Set::Scalar /*time*/, Set::Scalar dt)
98 : {
99 : //
100 : // FFT Boilerplate
101 : //
102 : amrex::FFT::R2C my_fft(this->geom[lev].Domain());
103 : auto const &[cba, cdm] = my_fft.getSpectralDataLayout();
104 : const Set::Scalar* DX = geom[lev].CellSize();
105 : amrex::Box const & domain = this->geom[lev].Domain();
106 : Set::Scalar
107 : AMREX_D_DECL(
108 : pi_Lx = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(0) / DX[0],
109 : pi_Ly = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(1) / DX[1],
110 : pi_Lz = 2.0 * Set::Constant::Pi / geom[lev].Domain().length(2) / DX[2]);
111 : Set::Scalar scaling = 1.0 / geom[lev].Domain().d_numPts();
112 :
113 :
114 : //
115 : // Compute the gradient of the chemical potential in realspace
116 : //
117 : for ( amrex::MFIter mfi(*etanew_mf[lev],true); mfi.isValid(); ++mfi )
118 : {
119 : const amrex::Box& bx = mfi.tilebox();
120 : amrex::Array4<const amrex::Real> const& eta = etanew_mf[lev]->array(mfi);
121 : amrex::Array4<amrex::Real> const& inter = intermediate[lev]->array(mfi);
122 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
123 : {
124 : inter(i,j,k) = eta(i,j,k)*eta(i,j,k)*eta(i,j,k) - eta(i,j,k);
125 : });
126 : }
127 :
128 : intermediate[lev]->FillBoundary();
129 :
130 : //
131 : // FFT of eta
132 : //
133 : amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > eta_hat_mf(cba, cdm, 1, 0);
134 : my_fft.forward(*etanew_mf[lev], eta_hat_mf);
135 :
136 : //
137 : // FFT of chemical potential gradient
138 : //
139 : amrex::FabArray<amrex::BaseFab<amrex::GpuComplex<Set::Scalar> > > chempot_hat_mf(cba, cdm, 1, 0);
140 : my_fft.forward(*intermediate[lev], chempot_hat_mf);
141 :
142 : //
143 : // Perform update in spectral coordinatees
144 : //
145 : //for (amrex::MFIter mfi(eta_hat_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
146 : for (amrex::MFIter mfi(eta_hat_mf, false); mfi.isValid(); ++mfi)
147 : {
148 : const amrex::Box &bx = mfi.tilebox();
149 :
150 : Util::Message(INFO, bx);
151 :
152 : amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & eta_hat = eta_hat_mf.array(mfi);
153 : amrex::Array4<amrex::GpuComplex<Set::Scalar>> const & chempot_hat = chempot_hat_mf.array(mfi);
154 :
155 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int m, int n, int p) {
156 :
157 :
158 : // Get spectral coordinates
159 : AMREX_D_TERM(
160 : Set::Scalar k1 = m * pi_Lx;,
161 : Set::Scalar k2 = (n < domain.length(1)/2 ? n * pi_Ly : (n - domain.length(1)) * pi_Ly);,
162 : Set::Scalar k3 = (p < domain.length(2)/2 ? p * pi_Lz : (p - domain.length(2)) * pi_Lz););
163 :
164 :
165 : Set::Scalar lap = AMREX_D_TERM(k1 * k1, + k2 * k2, + k3*k3);
166 :
167 : Set::Scalar bilap = lap*lap;
168 :
169 : eta_hat(m, n, p) =
170 : (eta_hat(m, n, p) - L * lap * chempot_hat(m, n, p) * dt) /
171 : (1.0 + L * gamma * bilap * dt);
172 :
173 :
174 : eta_hat(m,n,p) *= scaling;
175 : });
176 : }
177 :
178 : //
179 : // Transform solution back to realspace
180 : //
181 : my_fft.backward(eta_hat_mf, *etanew_mf[lev]);
182 : }
183 : #else
184 : void
185 0 : CahnHilliard::AdvanceSpectral (int, Set::Scalar, Set::Scalar)
186 : {
187 0 : Util::Abort(INFO,"Alamo must be compiled with fft");
188 0 : }
189 : #endif
190 :
191 :
192 :
193 : void
194 0 : CahnHilliard::Initialize (int lev)
195 : {
196 0 : intermediate[lev]->setVal(0.0);
197 0 : ic->Initialize(lev,etanew_mf);
198 0 : if (method == "realspace")
199 0 : ic->Initialize(lev,etaold_mf);
200 0 : }
201 :
202 :
203 : void
204 0 : CahnHilliard::TagCellsForRefinement (int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/)
205 : {
206 0 : if (method == "spectral") return;
207 :
208 0 : const Set::Scalar* DX = geom[lev].CellSize();
209 0 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
210 :
211 0 : for (amrex::MFIter mfi(*etanew_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
212 : {
213 0 : const amrex::Box& bx = mfi.tilebox();
214 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
215 0 : Set::Patch<const Set::Scalar> eta = (*etanew_mf[lev]).array(mfi);
216 :
217 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
218 : {
219 0 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
220 0 : if (grad.lpNorm<2>() * dr > refinement_threshold)
221 0 : tags(i, j, k) = amrex::TagBox::SET;
222 0 : });
223 0 : }
224 : }
225 :
226 :
227 : }
|
