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