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