Line data Source code
1 : // TODO: Remove these
2 :
3 : #include "Elastic.H"
4 : #include "Set/Set.H"
5 :
6 : #include "Numeric/Stencil.H"
7 : namespace Operator
8 : {
9 : template<int SYM>
10 734 : Elastic<SYM>::Elastic(const Vector<Geometry>& a_geom,
11 : const Vector<BoxArray>& a_grids,
12 : const Vector<DistributionMapping>& a_dmap,
13 734 : const LPInfo& a_info)
14 : {
15 : BL_PROFILE("Operator::Elastic::Elastic()");
16 :
17 734 : define(a_geom, a_grids, a_dmap, a_info);
18 734 : }
19 :
20 : template<int SYM>
21 734 : Elastic<SYM>::~Elastic()
22 734 : {}
23 :
24 : template<int SYM>
25 : void
26 734 : Elastic<SYM>::define(const Vector<Geometry>& a_geom,
27 : const Vector<BoxArray>& a_grids,
28 : const Vector<DistributionMapping>& a_dmap,
29 : const LPInfo& a_info,
30 : const Vector<FabFactory<FArrayBox> const*>& a_factory)
31 : {
32 : BL_PROFILE("Operator::Elastic::define()");
33 :
34 734 : Operator::define(a_geom, a_grids, a_dmap, a_info, a_factory);
35 :
36 734 : int model_nghost = 2;
37 :
38 734 : m_ddw_mf.resize(m_num_amr_levels);
39 734 : m_psi_mf.resize(m_num_amr_levels);
40 1522 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
41 : {
42 788 : m_ddw_mf[amrlev].resize(m_num_mg_levels[amrlev]);
43 788 : m_psi_mf[amrlev].resize(m_num_mg_levels[amrlev]);
44 2395 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
45 : {
46 4821 : m_ddw_mf[amrlev][mglev].reset(new MultiTab(amrex::convert(m_grids[amrlev][mglev],
47 3214 : amrex::IntVect::TheNodeVector()),
48 : m_dmap[amrlev][mglev], 1, model_nghost));
49 1607 : m_psi_mf[amrlev][mglev].reset(new MultiFab(m_grids[amrlev][mglev],
50 : m_dmap[amrlev][mglev], 1, model_nghost));
51 :
52 1607 : if (!m_psi_set) m_psi_mf[amrlev][mglev]->setVal(1.0);
53 : }
54 : }
55 734 : }
56 :
57 : template <int SYM>
58 : void
59 0 : Elastic<SYM>::SetModel(MATRIX4& a_model)
60 : {
61 0 : for (int amrlev = 0; amrlev < m_num_amr_levels; amrlev++)
62 : {
63 0 : amrex::Box domain(m_geom[amrlev][0].Domain());
64 0 : domain.convert(amrex::IntVect::TheNodeVector());
65 :
66 0 : for (MFIter mfi(*m_ddw_mf[amrlev][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
67 : {
68 0 : Box bx = mfi.grownnodaltilebox();
69 :
70 0 : amrex::Array4<MATRIX4> const& ddw = (*(m_ddw_mf[amrlev][0])).array(mfi);
71 :
72 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
73 0 : ddw(i, j, k) = a_model;
74 :
75 : #ifdef AMREX_DEBUG
76 : if (ddw(i, j, k).contains_nan()) Util::Abort(INFO, "model is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev);
77 : #endif
78 : });
79 : }
80 :
81 0 : (*(m_ddw_mf[amrlev][0])).setMultiGhost(true);
82 0 : (*(m_ddw_mf[amrlev][0])).FillBoundary( Geom(amrlev,0).periodicity());
83 : }
84 0 : m_model_set = true;
85 0 : }
86 :
87 : template <int SYM>
88 : void
89 1688 : Elastic<SYM>::SetModel(int amrlev, const amrex::FabArray<amrex::BaseFab<MATRIX4> >& a_model)
90 : {
91 : BL_PROFILE("Operator::Elastic::SetModel()");
92 :
93 1688 : amrex::Box domain(m_geom[amrlev][0].Domain());
94 1688 : domain.convert(amrex::IntVect::TheNodeVector());
95 :
96 1688 : if (a_model.boxArray() != m_ddw_mf[amrlev][0]->boxArray()) Util::Abort(INFO, "Inconsistent box arrays\n", "a_model.boxArray()=\n", a_model.boxArray(), "\n but the current box array is \n", m_ddw_mf[amrlev][0]->boxArray());
97 1688 : if (a_model.DistributionMap() != m_ddw_mf[amrlev][0]->DistributionMap()) Util::Abort(INFO, "Inconsistent distribution maps");
98 1688 : if (a_model.nComp() != m_ddw_mf[amrlev][0]->nComp()) Util::Abort(INFO, "Inconsistent # of components - should be ", m_ddw_mf[amrlev][0]->nComp());
99 1688 : if (a_model.nGrow() != m_ddw_mf[amrlev][0]->nGrow()) Util::Abort(INFO, "Inconsistent # of ghost nodes, should be ", m_ddw_mf[amrlev][0]->nGrow());
100 :
101 :
102 3509 : for (MFIter mfi(a_model, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
103 : {
104 1821 : Box bx = mfi.grownnodaltilebox();
105 :
106 1821 : amrex::Array4<MATRIX4> const& C = (*(m_ddw_mf[amrlev][0])).array(mfi);
107 1821 : amrex::Array4<const MATRIX4> const& a_C = a_model.array(mfi);
108 :
109 1406399 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
110 1404578 : C(i, j, k) = a_C(i, j, k);
111 : });
112 : }
113 1688 : m_ddw_mf[amrlev][0]->setMultiGhost(true);
114 1688 : m_ddw_mf[amrlev][0]->FillBoundaryAndSync(Geom(amrlev,0).periodicity());
115 1688 : m_model_set = true;
116 1688 : }
117 :
118 : template <int SYM>
119 : void
120 44 : Elastic<SYM>::SetPsi(int amrlev, const amrex::MultiFab& a_psi_mf)
121 : {
122 : BL_PROFILE("Operator::Elastic::SetPsi()");
123 44 : amrex::Box domain(m_geom[amrlev][0].Domain());
124 :
125 158 : for (MFIter mfi(a_psi_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
126 : {
127 114 : Box bx = mfi.growntilebox() & domain;
128 :
129 114 : amrex::Array4<Set::Scalar> const& m_psi = (*(m_psi_mf[amrlev][0])).array(mfi);
130 114 : amrex::Array4<const Set::Scalar> const& a_psi = a_psi_mf.array(mfi);
131 :
132 165298 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
133 165184 : m_psi(i, j, k) = a_psi(i, j, k);
134 : });
135 : }
136 44 : m_psi_set = true;
137 44 : }
138 :
139 : template<int SYM>
140 : void
141 247822 : Elastic<SYM>::Fapply(int amrlev, int mglev, MultiFab& a_f, const MultiFab& a_u) const
142 : {
143 : BL_PROFILE("Operator::Elastic::Fapply()");
144 :
145 247822 : amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
146 247822 : domain.convert(amrex::IntVect::TheNodeVector());
147 :
148 247822 : amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(2));
149 247822 : stencilbox.convert(amrex::IntVect::TheNodeVector());
150 :
151 247822 : const Real* DX = m_geom[amrlev][mglev].CellSize();
152 :
153 509875 : for (MFIter mfi(a_f, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
154 : {
155 262053 : Box bx = mfi.validbox().grow(1) & domain;
156 262053 : amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
157 :
158 262053 : amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
159 262053 : amrex::Array4<const amrex::Real> const& U = a_u.array(mfi);
160 262053 : amrex::Array4<amrex::Real> const& F = a_f.array(mfi);
161 262053 : amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
162 :
163 262053 : const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
164 :
165 44790421 : amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
166 :
167 22264184 : Set::Vector f = Set::Vector::Zero();
168 :
169 22264184 : Set::Vector u;
170 73868922 : for (int p = 0; p < AMREX_SPACEDIM; p++) u(p) = U(i, j, k, p);
171 :
172 :
173 22264184 : bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
174 22264184 : AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
175 :
176 : // Determine if a special stencil will be necessary for first derivatives
177 : std::array<Numeric::StencilType, AMREX_SPACEDIM>
178 22264184 : sten = Numeric::GetStencil(i, j, k, stencilbox);
179 :
180 : // The displacement gradient tensor
181 22264184 : Set::Matrix gradu; // gradu(i,j) = u_{i,j)
182 :
183 : // Fill gradu
184 73868922 : for (int p = 0; p < AMREX_SPACEDIM; p++)
185 : {
186 176043324 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(U, i, j, k, p, DX, sten));,
187 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(U, i, j, k, p, DX, sten));,
188 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(U, i, j, k, p, DX, sten)););
189 : }
190 :
191 22264184 : Set::Scalar psi_avg = 1.0;
192 26801203 : if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
193 :
194 : // Stress tensor computed using the model fab
195 66792552 : Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
196 :
197 : // Boundary conditions
198 : /// \todo Important: we need a way to handle corners and edges.
199 22264184 : amrex::IntVect m(AMREX_D_DECL(i, j, k));
200 22264184 : if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
201 : {
202 7793281 : f = (*m_bc)(u, gradu, sig, i, j, k, stencilbox);
203 : }
204 : else
205 : {
206 :
207 :
208 : // The gradient of the displacement gradient tensor
209 : // TODO - replace with this call. But not for this PR
210 : //Set::Matrix3 gradgradu = Numeric::Hessian(U,i,j,k,DX,sten); // gradgradu[k](l,j) = u_{k,lj}
211 14470903 : Set::Matrix3 gradgradu; // gradgradu[k](l,j) = u_{k,lj}
212 :
213 : // Fill gradu and gradgradu
214 44456507 : for (int p = 0; p < AMREX_SPACEDIM; p++)
215 : {
216 : // Diagonal terms:
217 126205204 : AMREX_D_TERM(
218 : gradgradu(p, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(U, i, j, k, p, DX));,
219 : gradgradu(p, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(U, i, j, k, p, DX));,
220 : gradgradu(p, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(U, i, j, k, p, DX)););
221 :
222 : // Off-diagonal terms:
223 174979172 : AMREX_D_TERM(
224 : ,// 2D
225 : gradgradu(p, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(U, i, j, k, p, DX));
226 : gradgradu(p, 1, 0) = gradgradu(p, 0, 1);
227 : ,// 3D
228 : gradgradu(p, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(U, i, j, k, p, DX));
229 : gradgradu(p, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(U, i, j, k, p, DX));
230 : gradgradu(p, 2, 0) = gradgradu(p, 0, 2);
231 : gradgradu(p, 2, 1) = gradgradu(p, 1, 2););
232 : }
233 :
234 : //
235 : // Operator
236 : //
237 : // The return value is
238 : // f = C(grad grad u) + grad(C)*grad(u)
239 : // In index notation
240 : // f_i = C_{ijkl,j} u_{k,l} + C_{ijkl}u_{k,lj}
241 : //
242 :
243 43412709 : f = (DDW(i, j, k) * gradgradu) * psi_avg;
244 :
245 14470903 : if (!m_uniform)
246 : {
247 : MATRIX4
248 43753202 : AMREX_D_DECL(Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(DDW, i, j, k, 0, DX, sten)),
249 : Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(DDW, i, j, k, 0, DX, sten)),
250 : Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(DDW, i, j, k, 0, DX, sten)));
251 30688909 : f += (AMREX_D_TERM((Cgrad1 * gradu).col(0),
252 : +(Cgrad2 * gradu).col(1),
253 : +(Cgrad3 * gradu).col(2))) * (psi_avg);
254 : }
255 14470903 : if (m_psi_set)
256 : {
257 4281099 : Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX);
258 4281099 : gradpsi *= (1.0 - m_psi_small);
259 12843297 : f += (DDW(i, j, k) * gradu) * gradpsi;
260 : }
261 :
262 14470903 : if (std::isnan(f(0)) || std::isnan(f(1)))
263 : {
264 0 : Util::Message(INFO," ================= ");
265 0 : Util::Message(INFO,"amrlev=",amrlev);
266 0 : Util::Message(INFO,"mglev=",mglev);
267 0 : Util::Message(INFO,"i=",i," j=",j);
268 0 : Util::Message(INFO,"f: ",f.transpose());
269 0 : Util::Message(INFO,"U(i,j,k): ",U(i,j,k,0)," ",U(i,j,k,1));
270 0 : Util::Message(INFO,"U(i-1,j,k): ",U(i-1,j,k,0)," ",U(i-1,j,k,1));
271 0 : Util::Message(INFO,"U(i+1,j,k): ",U(i+1,j,k,0)," ",U(i-1,j,k,1));
272 0 : Util::Message(INFO,"U(i,j-1,k): ",U(i,j-1,k,0)," ",U(i,j-1,k,1));
273 0 : Util::Message(INFO,"U(i,j+1,k): ",U(i,j+1,k,0)," ",U(i,j+1,k,1));
274 0 : Util::Message(INFO,"gradu: ",gradu);
275 0 : Util::Message(INFO,"gradgradu[0]: ",gradgradu[0]);
276 0 : Util::Message(INFO,"gradgradu[1]: ",gradgradu[1]);
277 0 : Util::Message(INFO,"DDW (i ,j ): ",DDW(i,j,k));
278 0 : Util::Message(INFO,"DDW (i-1,j ): ",DDW(i-1,j,k));
279 0 : Util::Message(INFO,"DDW (i+1,j ): ",DDW(i+1,j,k));
280 0 : Util::Message(INFO,"DDW (i ,j+1): ",DDW(i,j+1,k));
281 0 : Util::Message(INFO,"DDW (i ,j-1): ",DDW(i,j-1,k));
282 0 : Util::Message(INFO,"psi_av: ",psi_avg);
283 0 : Util::Message(INFO,"psi_set: ",m_psi_set);
284 0 : Util::Message(INFO," ================= ");
285 0 : Util::Abort(INFO);
286 : }
287 :
288 : }
289 73868922 : AMREX_D_TERM(F(i, j, k, 0) = f[0];, F(i, j, k, 1) = f[1];, F(i, j, k, 2) = f[2];);
290 : });
291 : }
292 247822 : }
293 :
294 :
295 :
296 : template<int SYM>
297 : void
298 1607 : Elastic<SYM>::Diagonal(int amrlev, int mglev, MultiFab& a_diag)
299 : {
300 : BL_PROFILE("Operator::Elastic::Diagonal()");
301 :
302 1607 : amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
303 1607 : domain.convert(amrex::IntVect::TheNodeVector());
304 :
305 1607 : amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(2));
306 1607 : stencilbox.convert(amrex::IntVect::TheNodeVector());
307 :
308 1607 : const Real* DX = m_geom[amrlev][mglev].CellSize();
309 :
310 3347 : for (MFIter mfi(a_diag, false); mfi.isValid(); ++mfi)
311 : {
312 1740 : Box bx = mfi.validbox().grow(1) & domain;
313 1740 : amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
314 :
315 1740 : amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
316 1740 : amrex::Array4<Set::Scalar> const& diag = a_diag.array(mfi);
317 1740 : amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
318 :
319 1740 : const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
320 :
321 264229 : amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
322 : {
323 :
324 : // Determine if a special stencil will be necessary for first derivatives
325 : std::array<Numeric::StencilType, AMREX_SPACEDIM>
326 262489 : sten = Numeric::GetStencil(i, j, k, stencilbox);
327 :
328 : // gradu(i,j) = u_{i,j)
329 262489 : std::array<Set::Matrix,AMREX_SPACEDIM> gradu = Numeric::Gradient_Diagonal<Set::Matrix>(DX, sten);
330 :
331 : // gradgradu[k](l,j) = u_{k,lj}
332 262489 : std::array<Set::Matrix3,AMREX_SPACEDIM> gradgradu = Numeric::Gradient_Diagonal<Set::Matrix3>(DX);
333 :
334 :
335 262489 : Set::Vector f = Set::Vector::Zero();
336 :
337 : bool
338 262489 : AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
339 262489 : AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
340 :
341 262489 : Set::Scalar psi_avg = 1.0;
342 344205 : if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
343 :
344 :
345 :
346 878667 : for (int p = 0; p < AMREX_SPACEDIM; p++)
347 : {
348 :
349 616178 : diag(i, j, k, p) = 0.0;
350 :
351 :
352 616178 : amrex::IntVect m(AMREX_D_DECL(i, j, k));
353 616178 : if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
354 : {
355 722610 : Set::Matrix sig = DDW(i, j, k) * gradu[p] * psi_avg;
356 240870 : Set::Vector u = Set::Vector::Zero();
357 240870 : u(p) = 1.0;
358 240870 : f = (*m_bc)(u, gradu[p], sig, i, j, k, stencilbox);
359 240870 : diag(i, j, k, p) = f(p);
360 240870 : }
361 : else
362 : {
363 1125924 : Set::Vector f = (DDW(i, j, k) * gradgradu[p]) * psi_avg;
364 750616 : diag(i, j, k, p) += f(p);
365 : }
366 :
367 : #ifdef AMREX_DEBUG
368 : if (std::isnan(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
369 : if (std::isinf(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is inf at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
370 : if (diag(i, j, k, p) == 0) Util::Abort(INFO, "diagonal is zero at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
371 : #endif
372 :
373 : }
374 : });
375 : }
376 :
377 1607 : a_diag.FillBoundaryAndSync(Geom(amrlev,mglev).periodicity());
378 1607 : nodalSync(amrlev,mglev,a_diag);
379 1607 : }
380 :
381 :
382 : template<int SYM>
383 : void
384 0 : Elastic<SYM>::Error0x(int amrlev, int mglev, MultiFab& R0x, const MultiFab& x) const
385 : {
386 : BL_PROFILE("Operator::Elastic::Error0x()");
387 0 : Util::Message(INFO);
388 :
389 0 : int ncomp = x.nComp();//getNComp();
390 0 : int nghost = x.nGrow();
391 :
392 0 : if (!m_diagonal_computed)
393 0 : Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
394 :
395 0 : amrex::MultiFab D0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
396 0 : amrex::MultiFab AD0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
397 :
398 0 : amrex::MultiFab::Copy(D0x, x, 0, 0, ncomp, nghost); // D0x = x
399 0 : amrex::MultiFab::Divide(D0x, *m_diag[amrlev][mglev], 0, 0, ncomp, 0); // D0x = x/diag
400 0 : amrex::MultiFab::Copy(AD0x, D0x, 0, 0, ncomp, nghost); // AD0x = D0x
401 :
402 0 : Fapply(amrlev, mglev, AD0x, D0x); // AD0x = A * D0 * x
403 :
404 0 : amrex::MultiFab::Copy(R0x, x, 0, 0, ncomp, nghost); // R0x = x
405 0 : amrex::MultiFab::Subtract(R0x, AD0x, 0, 0, ncomp, nghost); // R0x = x - AD0x
406 0 : }
407 :
408 :
409 : template<int SYM>
410 : void
411 0 : Elastic<SYM>::FFlux(int /*amrlev*/, const MFIter& /*mfi*/,
412 : const std::array<FArrayBox*, AMREX_SPACEDIM>& sigmafab,
413 : const FArrayBox& /*ufab*/, const int /*face_only*/) const
414 : {
415 : BL_PROFILE("Operator::Elastic::FFlux()");
416 0 : Util::Message(INFO);
417 0 : amrex::BaseFab<amrex::Real> AMREX_D_DECL(&fxfab = *sigmafab[0],
418 : &fyfab = *sigmafab[1],
419 : &fzfab = *sigmafab[2]);
420 0 : AMREX_D_TERM(fxfab.setVal(0.0);,
421 : fyfab.setVal(0.0);,
422 : fzfab.setVal(0.0););
423 :
424 0 : }
425 :
426 : template<int SYM>
427 : void
428 0 : Elastic<SYM>::Strain(int amrlev,
429 : amrex::MultiFab& a_eps,
430 : const amrex::MultiFab& a_u,
431 : bool voigt) const
432 : {
433 : BL_PROFILE("Operator::Elastic::Strain()");
434 :
435 0 : const amrex::Real* DX = m_geom[amrlev][0].CellSize();
436 0 : amrex::Box domain(m_geom[amrlev][0].Domain());
437 0 : domain.convert(amrex::IntVect::TheNodeVector());
438 :
439 :
440 0 : for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
441 : {
442 0 : const Box& bx = mfi.tilebox();
443 0 : amrex::Array4<amrex::Real> const& epsilon = a_eps.array(mfi);
444 0 : amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
445 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
446 : {
447 0 : Set::Matrix gradu;
448 :
449 : std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
450 0 : = Numeric::GetStencil(i, j, k, domain);
451 :
452 : // Fill gradu
453 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
454 : {
455 0 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
456 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
457 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
458 : }
459 :
460 0 : Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
461 :
462 0 : if (voigt)
463 : {
464 0 : AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
465 : ,
466 : epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(0, 1);
467 : ,
468 : epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(2, 2);
469 : epsilon(i, j, k, 3) = eps(1, 2); epsilon(i, j, k, 4) = eps(2, 0); epsilon(i, j, k, 5) = eps(0, 1););
470 : }
471 : else
472 : {
473 0 : AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
474 : ,
475 : epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1);
476 : epsilon(i, j, k, 2) = eps(1, 0); epsilon(i, j, k, 3) = eps(1, 1);
477 : ,
478 : epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1); epsilon(i, j, k, 2) = eps(0, 2);
479 : epsilon(i, j, k, 3) = eps(1, 0); epsilon(i, j, k, 4) = eps(1, 1); epsilon(i, j, k, 5) = eps(1, 2);
480 : epsilon(i, j, k, 6) = eps(2, 0); epsilon(i, j, k, 7) = eps(2, 1); epsilon(i, j, k, 8) = eps(2, 2););
481 : }
482 : });
483 : }
484 0 : }
485 :
486 :
487 : template<int SYM>
488 : void
489 0 : Elastic<SYM>::Stress(int amrlev,
490 : amrex::MultiFab& a_sigma,
491 : const amrex::MultiFab& a_u,
492 : bool voigt, bool a_homogeneous)
493 : {
494 : BL_PROFILE("Operator::Elastic::Stress()");
495 0 : SetHomogeneous(a_homogeneous);
496 :
497 0 : const amrex::Real* DX = m_geom[amrlev][0].CellSize();
498 0 : amrex::Box domain(m_geom[amrlev][0].Domain());
499 0 : domain.convert(amrex::IntVect::TheNodeVector());
500 :
501 0 : for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
502 : {
503 0 : const Box& bx = mfi.tilebox();
504 0 : amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
505 0 : amrex::Array4<amrex::Real> const& sigma = a_sigma.array(mfi);
506 0 : amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][0]->array(mfi);
507 0 : amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
508 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
509 : {
510 0 : Set::Matrix gradu;
511 :
512 : std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
513 0 : = Numeric::GetStencil(i, j, k, domain);
514 :
515 : // Fill gradu
516 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
517 : {
518 0 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
519 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
520 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
521 : }
522 :
523 0 : Set::Scalar psi_avg = 1.0;
524 0 : if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
525 0 : Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
526 :
527 0 : if (voigt)
528 : {
529 0 : AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
530 : ,
531 : sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(0, 1);
532 : ,
533 : sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(2, 2);
534 : sigma(i, j, k, 3) = sig(1, 2); sigma(i, j, k, 4) = sig(2, 0); sigma(i, j, k, 5) = sig(0, 1););
535 : }
536 : else
537 : {
538 0 : AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
539 : ,
540 : sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1);
541 : sigma(i, j, k, 2) = sig(1, 0); sigma(i, j, k, 3) = sig(1, 1);
542 : ,
543 : sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1); sigma(i, j, k, 2) = sig(0, 2);
544 : sigma(i, j, k, 3) = sig(1, 0); sigma(i, j, k, 4) = sig(1, 1); sigma(i, j, k, 5) = sig(1, 2);
545 : sigma(i, j, k, 6) = sig(2, 0); sigma(i, j, k, 7) = sig(2, 1); sigma(i, j, k, 8) = sig(2, 2););
546 : }
547 : });
548 : }
549 0 : }
550 :
551 :
552 : template<int SYM>
553 : void
554 0 : Elastic<SYM>::Energy(int amrlev,
555 : amrex::MultiFab& a_energy,
556 : const amrex::MultiFab& a_u, bool a_homogeneous)
557 : {
558 : BL_PROFILE("Operator::Elastic::Energy()");
559 0 : SetHomogeneous(a_homogeneous);
560 :
561 0 : amrex::Box domain(m_geom[amrlev][0].Domain());
562 0 : domain.convert(amrex::IntVect::TheNodeVector());
563 :
564 0 : const amrex::Real* DX = m_geom[amrlev][0].CellSize();
565 :
566 0 : for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
567 : {
568 0 : const Box& bx = mfi.tilebox();
569 0 : amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
570 0 : amrex::Array4<amrex::Real> const& energy = a_energy.array(mfi);
571 0 : amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
572 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
573 : {
574 0 : Set::Matrix gradu;
575 :
576 : std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
577 0 : = Numeric::GetStencil(i, j, k, domain);
578 :
579 : // Fill gradu
580 0 : for (int p = 0; p < AMREX_SPACEDIM; p++)
581 : {
582 0 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
583 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
584 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
585 : }
586 :
587 0 : Set::Matrix eps = .5 * (gradu + gradu.transpose());
588 0 : Set::Matrix sig = DDW(i, j, k) * gradu;
589 :
590 : // energy(i,j,k) = (gradu.transpose() * sig).trace();
591 :
592 : //Util::Abort(INFO,"Fix this"); //
593 : //energy(i,j,k) = C(i,j,k).W(gradu);
594 0 : for (int m = 0; m < AMREX_SPACEDIM; m++)
595 : {
596 0 : for (int n = 0; n < AMREX_SPACEDIM; n++)
597 : {
598 0 : energy(i, j, k) += .5 * sig(m, n) * eps(m, n);
599 : }
600 : }
601 : });
602 : }
603 0 : }
604 :
605 : template<int SYM>
606 : void
607 734 : Elastic<SYM>::averageDownCoeffs()
608 : {
609 : BL_PROFILE("Elastic::averageDownCoeffs()");
610 :
611 734 : if (m_average_down_coeffs)
612 0 : for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
613 0 : averageDownCoeffsDifferentAmrLevels(amrlev);
614 :
615 734 : averageDownCoeffsSameAmrLevel(0);
616 1522 : for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
617 : {
618 2395 : for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
619 : {
620 1607 : if (m_ddw_mf[amrlev][mglev]) {
621 1607 : FillBoundaryCoeff(*m_ddw_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
622 1607 : FillBoundaryCoeff(*m_psi_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
623 : }
624 : }
625 : }
626 734 : }
627 :
628 : template<int SYM>
629 : void
630 0 : Elastic<SYM>::averageDownCoeffsDifferentAmrLevels(int fine_amrlev)
631 : {
632 : BL_PROFILE("Operator::Elastic::averageDownCoeffsDifferentAmrLevels()");
633 0 : Util::Assert(INFO, TEST(fine_amrlev > 0));
634 :
635 0 : const int crse_amrlev = fine_amrlev - 1;
636 0 : const int ncomp = 1;
637 :
638 0 : MultiTab& crse_ddw = *m_ddw_mf[crse_amrlev][0];
639 0 : MultiTab& fine_ddw = *m_ddw_mf[fine_amrlev][0];
640 :
641 0 : amrex::Box cdomain(m_geom[crse_amrlev][0].Domain());
642 0 : cdomain.convert(amrex::IntVect::TheNodeVector());
643 :
644 0 : const Geometry& cgeom = m_geom[crse_amrlev][0];
645 :
646 0 : const BoxArray& fba = fine_ddw.boxArray();
647 0 : const DistributionMapping& fdm = fine_ddw.DistributionMap();
648 :
649 0 : MultiTab fine_ddw_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
650 0 : fine_ddw_for_coarse.ParallelCopy(crse_ddw, 0, 0, ncomp, 0, 0, cgeom.periodicity());
651 :
652 0 : const int coarse_fine_node = 1;
653 0 : const int fine_fine_node = 2;
654 :
655 0 : amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
656 0 : nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
657 :
658 0 : amrex::iMultiFab cellmask(amrex::convert(amrex::coarsen(fba, 2), amrex::IntVect::TheCellVector()), fdm, 1, 2);
659 0 : cellmask.ParallelCopy(*m_cc_fine_mask[crse_amrlev], 0, 0, 1, 1, 1, cgeom.periodicity());
660 :
661 0 : for (MFIter mfi(fine_ddw_for_coarse, false); mfi.isValid(); ++mfi)
662 : {
663 0 : const Box& bx = mfi.validbox();
664 :
665 0 : amrex::Array4<const int> const& nmask = nodemask.array(mfi);
666 : //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
667 :
668 0 : amrex::Array4<MATRIX4> const& cdata = fine_ddw_for_coarse.array(mfi);
669 0 : amrex::Array4<const MATRIX4> const& fdata = fine_ddw.array(mfi);
670 :
671 0 : const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
672 :
673 0 : for (int n = 0; n < fine_ddw.nComp(); n++)
674 : {
675 : // I,J,K == coarse coordinates
676 : // i,j,k == fine coordinates
677 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
678 0 : int i = I * 2, j = J * 2, k = K * 2;
679 :
680 0 : if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
681 : {
682 0 : if ((I == lo.x || I == hi.x) &&
683 0 : (J == lo.y || J == hi.y) &&
684 0 : (K == lo.z || K == hi.z)) // Corner
685 0 : cdata(I, J, K, n) = fdata(i, j, k, n);
686 0 : else if ((J == lo.y || J == hi.y) &&
687 0 : (K == lo.z || K == hi.z)) // X edge
688 0 : cdata(I, J, K, n) = fdata(i - 1, j, k, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i + 1, j, k, n) * 0.25;
689 0 : else if ((K == lo.z || K == hi.z) &&
690 0 : (I == lo.x || I == hi.x)) // Y edge
691 0 : cdata(I, J, K, n) = fdata(i, j - 1, k, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i, j + 1, k, n) * 0.25;
692 0 : else if ((I == lo.x || I == hi.x) &&
693 0 : (J == lo.y || J == hi.y)) // Z edge
694 0 : cdata(I, J, K, n) = fdata(i, j, k - 1, n) * 0.25 + fdata(i, j, k, n) * 0.5 + fdata(i, j, k + 1, n) * 0.25;
695 0 : else if (I == lo.x || I == hi.x) // X face
696 0 : cdata(I, J, K, n) =
697 0 : (fdata(i, j - 1, k - 1, n) + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j + 1, k - 1, n)
698 0 : + fdata(i, j - 1, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j + 1, k, n) * 2.0
699 0 : + fdata(i, j - 1, k + 1, n) + fdata(i, j, k + 1, n) * 2.0 + fdata(i, j + 1, k + 1, n)) / 16.0;
700 0 : else if (J == lo.y || J == hi.y) // Y face
701 0 : cdata(I, J, K, n) =
702 0 : (fdata(i - 1, j, k - 1, n) + fdata(i - 1, j, k, n) * 2.0 + fdata(i - 1, j, k + 1, n)
703 0 : + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j, k + 1, n) * 2.0
704 0 : + fdata(i + 1, j, k - 1, n) + fdata(i + 1, j, k, n) * 2.0 + fdata(i + 1, j, k + 1, n)) / 16.0;
705 0 : else if (K == lo.z || K == hi.z) // Z face
706 0 : cdata(I, J, K, n) =
707 0 : (fdata(i - 1, j - 1, k, n) + fdata(i, j - 1, k, n) * 2.0 + fdata(i + 1, j - 1, k, n)
708 0 : + fdata(i - 1, j, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i + 1, j, k, n) * 2.0
709 0 : + fdata(i - 1, j + 1, k, n) + fdata(i, j + 1, k, n) * 2.0 + fdata(i + 1, j + 1, k, n)) / 16.0;
710 : else // Interior
711 0 : cdata(I, J, K, n) =
712 0 : (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
713 0 : fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
714 0 : +
715 0 : (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
716 0 : fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
717 0 : fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
718 0 : +
719 0 : (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
720 0 : fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
721 0 : +
722 0 : fdata(i, j, k, n) / 8.0;
723 :
724 : #ifdef AMREX_DEBUG
725 : if (cdata(I, J, K).contains_nan()) Util::Abort(INFO, "restricted model is nan at (", i, ",", j, ",", k, "), fine_amrlev=", fine_amrlev);
726 : #endif
727 : }
728 :
729 : });
730 : }
731 : }
732 :
733 : // Copy the fine residual restricted onto the coarse grid
734 : // into the final residual.
735 :
736 0 : crse_ddw.ParallelCopy(fine_ddw_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
737 : //const int mglev = 0;
738 : //Util::RealFillBoundary(crse_ddw, m_geom[crse_amrlev][mglev]);
739 :
740 0 : FillBoundaryCoeff(crse_ddw,Geom(fine_amrlev,0).periodicity());
741 :
742 0 : }
743 :
744 :
745 :
746 : template<int SYM>
747 : void
748 734 : Elastic<SYM>::averageDownCoeffsSameAmrLevel(int amrlev)
749 : {
750 : BL_PROFILE("Elastic::averageDownCoeffsSameAmrLevel()");
751 :
752 4691 : for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev)
753 : {
754 819 : amrex::Box cdomain(m_geom[amrlev][mglev].growPeriodicDomain(2));
755 819 : cdomain.convert(amrex::IntVect::TheNodeVector());
756 819 : amrex::Box fdomain(m_geom[amrlev][mglev - 1].Domain());
757 819 : fdomain.convert(amrex::IntVect::TheNodeVector());
758 :
759 819 : MultiTab& crse = *m_ddw_mf[amrlev][mglev];
760 819 : MultiTab& fine = *m_ddw_mf[amrlev][mglev - 1];
761 :
762 819 : amrex::BoxArray crseba = crse.boxArray();
763 819 : amrex::BoxArray fineba = fine.boxArray();
764 :
765 819 : BoxArray newba = crseba;
766 819 : newba.refine(2);
767 819 : MultiTab fine_on_crseba;
768 819 : fine_on_crseba.define(newba, crse.DistributionMap(), 1, 4);
769 819 : fine_on_crseba.ParallelCopy(fine, 0, 0, 1, 2, 4, m_geom[amrlev][mglev-1].periodicity());
770 : /* ine_on_crseba.FillBoundaryAndSync(m_geom[amrlev][mglev-1].periodicity()); */
771 :
772 1638 : for (MFIter mfi(crse, false); mfi.isValid(); ++mfi)
773 : {
774 :
775 819 : Box bx = mfi.grownnodaltilebox() & cdomain;
776 : /*Box bx = mfi.grownnodaltilebox(-1,1) & cdomain;*/
777 :
778 819 : amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>> const& fdata = fine_on_crseba.array(mfi);
779 819 : amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& cdata = crse.array(mfi);
780 :
781 819 : const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
782 : /*const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);*/
783 :
784 : // I,J,K == coarse coordinates
785 : // i,j,k == fine coordinates
786 57553 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
787 28367 : int i = 2 * I, j = 2 * J, k = 2 * K;
788 :
789 28367 : if ((I == lo.x || I == hi.x) &&
790 13426 : (J == lo.y || J == hi.y) &&
791 8076 : (K == lo.z || K == hi.z)) // Corner
792 17028 : cdata(I, J, K) = fdata(i, j, k);
793 22691 : else if ((J == lo.y || J == hi.y) &&
794 7646 : (K == lo.z || K == hi.z)) // X edge
795 38744 : cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
796 18645 : else if ((K == lo.z || K == hi.z) &&
797 13245 : (I == lo.x || I == hi.x)) // Y edge
798 39528 : cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
799 14495 : else if ((I == lo.x || I == hi.x) &&
800 3600 : (J == lo.y || J == hi.y)) // Z edge
801 24000 : cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
802 12095 : else if (I == lo.x || I == hi.x) // X face
803 1200 : cdata(I, J, K) =
804 4800 : (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
805 6000 : + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
806 18000 : + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
807 10895 : else if (J == lo.y || J == hi.y) // Y face
808 1200 : cdata(I, J, K) =
809 4800 : (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
810 6000 : + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
811 18000 : + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
812 9695 : else if (K == lo.z || K == hi.z) // Z face
813 9095 : cdata(I, J, K) =
814 48733 : (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
815 63307 : + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
816 78744 : + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
817 : else // Interior
818 600 : cdata(I, J, K) =
819 3400 : (fdata(i - 1, j - 1, k - 1) + fdata(i - 1, j - 1, k + 1) + fdata(i - 1, j + 1, k - 1) + fdata(i - 1, j + 1, k + 1) +
820 3000 : fdata(i + 1, j - 1, k - 1) + fdata(i + 1, j - 1, k + 1) + fdata(i + 1, j + 1, k - 1) + fdata(i + 1, j + 1, k + 1)) / 64.0
821 200 : +
822 3400 : (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
823 3200 : fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
824 3000 : fdata(i - 1, j - 1, k) + fdata(i - 1, j + 1, k) + fdata(i + 1, j - 1, k) + fdata(i + 1, j + 1, k)) / 32.0
825 200 : +
826 2600 : (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
827 2200 : fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
828 12000 : +
829 1600 : fdata(i, j, k) / 8.0;
830 :
831 : #ifdef AMREX_DEBUG
832 : if (cdata(I, J, K).contains_nan()) Util::Abort(INFO, "restricted model is nan at crse coordinates (I=", I, ",J=", J, ",K=", k, "), amrlev=", amrlev, " interpolating from mglev", mglev - 1, " to ", mglev);
833 : #endif
834 : });
835 : }
836 819 : FillBoundaryCoeff(crse, Geom(amrlev,mglev).periodicity());
837 :
838 :
839 819 : if (!m_psi_set) continue;
840 :
841 46 : amrex::Box cdomain_cell(m_geom[amrlev][mglev].Domain());
842 46 : amrex::Box fdomain_cell(m_geom[amrlev][mglev - 1].Domain());
843 46 : MultiFab& crse_psi = *m_psi_mf[amrlev][mglev];
844 46 : MultiFab& fine_psi = *m_psi_mf[amrlev][mglev - 1];
845 46 : MultiFab fine_psi_on_crseba;
846 92 : fine_psi_on_crseba.define(newba.convert(amrex::IntVect::TheCellVector()), crse_psi.DistributionMap(), 1, 1);
847 46 : fine_psi_on_crseba.ParallelCopy(fine_psi, 0, 0, 1, 1, 1, m_geom[amrlev][mglev].periodicity());
848 :
849 92 : for (MFIter mfi(crse_psi, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
850 : {
851 46 : Box bx = mfi.tilebox();
852 46 : bx = bx & cdomain_cell;
853 :
854 46 : amrex::Array4<const Set::Scalar> const& fdata = fine_psi_on_crseba.array(mfi);
855 46 : amrex::Array4<Set::Scalar> const& cdata = crse_psi.array(mfi);
856 :
857 46 : const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
858 :
859 : // I,J,K == coarse coordinates
860 : // i,j,k == fine coordinates
861 5806 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
862 2880 : int i = 2 * I, j = 2 * J, k = 2 * K;
863 :
864 2880 : if ((I == lo.x || I == hi.x) &&
865 260 : (J == lo.y || J == hi.y) &&
866 46 : (K == lo.z || K == hi.z)) // Corner
867 138 : cdata(I, J, K) = fdata(i, j, k);
868 2834 : else if ((J == lo.y || J == hi.y) &&
869 274 : (K == lo.z || K == hi.z)) // X edge
870 1370 : cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
871 2560 : else if ((K == lo.z || K == hi.z) &&
872 2560 : (I == lo.x || I == hi.x)) // Y edge
873 1070 : cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
874 2346 : else if ((I == lo.x || I == hi.x) &&
875 0 : (J == lo.y || J == hi.y)) // Z edge
876 0 : cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
877 2346 : else if (I == lo.x || I == hi.x) // X face
878 0 : cdata(I, J, K) =
879 0 : (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
880 0 : + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
881 0 : + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
882 2346 : else if (J == lo.y || J == hi.y) // Y face
883 0 : cdata(I, J, K) =
884 0 : (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
885 0 : + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
886 0 : + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
887 2346 : else if (K == lo.z || K == hi.z) // Z face
888 2346 : cdata(I, J, K) =
889 7038 : (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
890 7038 : + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
891 9384 : + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
892 : else // Interior
893 0 : cdata(I, J, K) =
894 0 : (fdata(i - 1, j - 1, k - 1) + fdata(i - 1, j - 1, k + 1) + fdata(i - 1, j + 1, k - 1) + fdata(i - 1, j + 1, k + 1) +
895 0 : fdata(i + 1, j - 1, k - 1) + fdata(i + 1, j - 1, k + 1) + fdata(i + 1, j + 1, k - 1) + fdata(i + 1, j + 1, k + 1)) / 64.0
896 0 : +
897 0 : (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
898 0 : fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
899 0 : fdata(i - 1, j - 1, k) + fdata(i - 1, j + 1, k) + fdata(i + 1, j - 1, k) + fdata(i + 1, j + 1, k)) / 32.0
900 0 : +
901 0 : (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
902 0 : fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
903 0 : +
904 0 : fdata(i, j, k) / 8.0;
905 : });
906 : }
907 46 : FillBoundaryCoeff(crse_psi, Geom(amrlev,mglev).periodicity());
908 :
909 : }
910 734 : }
911 :
912 : template<int SYM>
913 : void
914 2426 : Elastic<SYM>::FillBoundaryCoeff(MultiTab& sigma, const amrex::Periodicity& p)
915 : {
916 2426 : sigma.setMultiGhost(true);
917 2426 : sigma.FillBoundaryAndSync(p);
918 2426 : }
919 :
920 : template<int SYM>
921 : void
922 1653 : Elastic<SYM>::FillBoundaryCoeff(MultiFab& psi, const amrex::Periodicity& p)
923 : {
924 1653 : psi.setMultiGhost(true);
925 1653 : psi.FillBoundaryAndSync(p);
926 1653 : }
927 :
928 : template class Elastic<Set::Sym::Major>;
929 : template class Elastic<Set::Sym::Isotropic>;
930 : template class Elastic<Set::Sym::MajorMinor>;
931 : template class Elastic<Set::Sym::Diagonal>;
932 :
933 : }
|
