Alamo
Elastic.cpp
Go to the documentation of this file.
1// TODO: Remove these
2
3#include "Elastic.H"
4#include "Set/Set.H"
5
6#include "Numeric/Stencil.H"
7namespace Operator
8{
9template<int SYM>
10Elastic<SYM>::Elastic(const Vector<Geometry>& a_geom,
11 const Vector<BoxArray>& a_grids,
12 const Vector<DistributionMapping>& a_dmap,
13 const LPInfo& a_info,
14 bool a_conservative_face_flux)
15{
16 BL_PROFILE("Operator::Elastic::Elastic()");
17
18 define(a_geom, a_grids, a_dmap, a_info, {},
19 a_conservative_face_flux);
20}
21
22template<int SYM>
25
26template<int SYM>
27void
28Elastic<SYM>::define(const Vector<Geometry>& a_geom,
29 const Vector<BoxArray>& a_grids,
30 const Vector<DistributionMapping>& a_dmap,
31 const LPInfo& a_info,
32 const Vector<FabFactory<FArrayBox> const*>& a_factory,
33 bool a_conservative_face_flux)
34{
35 BL_PROFILE("Operator::Elastic::define()");
36
37 m_conservative_face_flux = a_conservative_face_flux;
38 Operator::define(a_geom, a_grids, a_dmap, a_info, a_factory);
39
40 int model_nghost = 2;
41 // A cell-to-node average at the outer diagonal ghost row reaches one
42 // cell farther than the nodal coefficient stencil.
43 int psi_nghost = model_nghost + 1;
44 int model_ncomp = m_conservative_face_flux ? AMREX_SPACEDIM + 1 : 1;
45
46 m_ddw_mf.resize(m_num_amr_levels);
47 m_psi_mf.resize(m_num_amr_levels);
48 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
49 {
50 m_ddw_mf[amrlev].resize(m_num_mg_levels[amrlev]);
51 m_psi_mf[amrlev].resize(m_num_mg_levels[amrlev]);
52 for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
53 {
54 m_ddw_mf[amrlev][mglev].reset(new MultiTab(amrex::convert(m_grids[amrlev][mglev],
55 amrex::IntVect::TheNodeVector()),
56 m_dmap[amrlev][mglev], model_ncomp, model_nghost));
57 m_psi_mf[amrlev][mglev].reset(new MultiFab(m_grids[amrlev][mglev],
58 m_dmap[amrlev][mglev], 1, psi_nghost));
59
60 if (!m_psi_set) m_psi_mf[amrlev][mglev]->setVal(1.0);
61 }
62 }
63}
64
65template <int SYM>
66void
67Elastic<SYM>::SetModel(MATRIX4& a_model)
68{
69 for (int amrlev = 0; amrlev < m_num_amr_levels; amrlev++)
70 {
71 amrex::Box domain(m_geom[amrlev][0].Domain());
72 domain.convert(amrex::IntVect::TheNodeVector());
73
74 for (MFIter mfi(*m_ddw_mf[amrlev][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
75 {
76 Box bx = mfi.grownnodaltilebox();
77
78 amrex::Array4<MATRIX4> const& ddw = (*(m_ddw_mf[amrlev][0])).array(mfi);
79 const int ncomp = m_ddw_mf[amrlev][0]->nComp();
80
81 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
82 for (int n = 0; n < ncomp; ++n)
83 ddw(i, j, k, n) = a_model;
84
85#ifdef AMREX_DEBUG
86 if (ddw(i, j, k, 0).contains_nan()) Util::Abort(INFO, "model is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev);
87#endif
88 });
89 }
90
91 (*(m_ddw_mf[amrlev][0])).setMultiGhost(true);
92 (*(m_ddw_mf[amrlev][0])).FillBoundary( Geom(amrlev,0).periodicity());
93 }
94 m_model_set = true;
95}
96
97template <int SYM>
98void
99Elastic<SYM>::SetModel(int amrlev, const amrex::FabArray<amrex::BaseFab<MATRIX4> >& a_model)
100{
101 BL_PROFILE("Operator::Elastic::SetModel()");
102
103 amrex::Box domain(m_geom[amrlev][0].Domain());
104 domain.convert(amrex::IntVect::TheNodeVector());
105
106 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());
107 if (a_model.DistributionMap() != m_ddw_mf[amrlev][0]->DistributionMap()) Util::Abort(INFO, "Inconsistent distribution maps");
108 if (a_model.nComp() != 1 &&
109 a_model.nComp() != m_ddw_mf[amrlev][0]->nComp())
110 Util::Abort(INFO, "Inconsistent # of coefficient components - should be 1 or ",
111 m_ddw_mf[amrlev][0]->nComp());
112 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());
113
114 const bool nodal_only = a_model.nComp() == 1;
115 const int ncomp = m_ddw_mf[amrlev][0]->nComp();
116
117 for (MFIter mfi(a_model, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
118 {
119 Box bx = mfi.grownnodaltilebox();
120
121 amrex::Array4<MATRIX4> const& C = (*(m_ddw_mf[amrlev][0])).array(mfi);
122 amrex::Array4<const MATRIX4> const& a_C = a_model.array(mfi);
123
124 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
125 C(i, j, k, 0) = a_C(i, j, k, 0);
126 for (int n = 1; n < ncomp; ++n)
127 C(i, j, k, n) = a_C(i, j, k, nodal_only ? 0 : n);
128 });
129 }
130 m_ddw_mf[amrlev][0]->setMultiGhost(true);
131 m_ddw_mf[amrlev][0]->FillBoundaryAndSync(Geom(amrlev,0).periodicity());
132 m_model_set = true;
133}
134
135template <int SYM>
136void
137Elastic<SYM>::SetPsi(int amrlev, const amrex::MultiFab& a_psi_mf)
138{
139 BL_PROFILE("Operator::Elastic::SetPsi()");
140 amrex::Box domain(m_geom[amrlev][0].Domain());
141
142 for (MFIter mfi(a_psi_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
143 {
144 Box bx = mfi.growntilebox() & domain;
145
146 amrex::Array4<Set::Scalar> const& m_psi = (*(m_psi_mf[amrlev][0])).array(mfi);
147 amrex::Array4<const Set::Scalar> const& a_psi = a_psi_mf.array(mfi);
148
149 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
150 m_psi(i, j, k) = a_psi(i, j, k);
151 });
152 }
153 m_psi_set = true;
154}
155
156template<int SYM>
157void
158Elastic<SYM>::Fapply(int amrlev, int mglev, MultiFab& a_f, const MultiFab& a_u) const
159{
160 BL_PROFILE("Operator::Elastic::Fapply()");
161
162 amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
163 domain.convert(amrex::IntVect::TheNodeVector());
164
165 amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(2));
166 stencilbox.convert(amrex::IntVect::TheNodeVector());
167
168 const Real* DX = m_geom[amrlev][mglev].CellSize();
169
170 for (MFIter mfi(a_f, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
171 {
172 Box bx = mfi.validbox().grow(1) & domain;
173 amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
174
175 amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
176 amrex::Array4<const amrex::Real> const& U = a_u.array(mfi);
177 amrex::Array4<amrex::Real> const& F = a_f.array(mfi);
178 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
179
180 if (m_conservative_face_flux)
181 {
182 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
183 amrex::LoopConcurrentOnCpu(tilebox, [=] (int i, int j, int k)
184 {
185 Set::Vector f = Set::Vector::Zero();
186 Set::Vector u;
187 for (int p = 0; p < AMREX_SPACEDIM; ++p)
188 u(p) = U(i, j, k, p);
189
190 const auto sten = Numeric::GetStencil(i, j, k, stencilbox);
191 Set::Matrix gradu = Numeric::Gradient(U, i, j, k, DX, sten);
192 const int index[3] = {i, j, k};
193 const int lower[3] = {lo.x, lo.y, lo.z};
194 const int upper[3] = {hi.x, hi.y, hi.z};
195 bool on_boundary = false;
196 for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
197 on_boundary = on_boundary ||
198 index[dir] == lower[dir] || index[dir] == upper[dir];
199
200 if (on_boundary)
201 {
202 Set::Scalar psi_avg = 1.0;
203 if (m_psi_set)
204 psi_avg = (1.0 - m_psi_small) *
205 Numeric::Interpolate::CellToNodeAverage(
206 psi, i, j, k, 0) + m_psi_small;
207 const Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
208 f = (*m_bc)(u, gradu, sig, i, j, k, stencilbox);
209 }
210 else
211 {
212 for (int face = 0; face < AMREX_SPACEDIM; ++face)
213 {
214 const int im = i - (face == 0);
215 const int jm = j - (face == 1);
216 const int km = k - (face == 2);
217 const Set::Matrix grad_hi =
218 Numeric::FaceGradient(U, i, j, k, face, DX);
219 const Set::Matrix grad_lo =
220 Numeric::FaceGradient(U, im, jm, km, face, DX);
221 const Set::Matrix flux_hi =
222 DDW(i, j, k, face + 1) * grad_hi;
223 const Set::Matrix flux_lo =
224 DDW(im, jm, km, face + 1) * grad_lo;
225 f += (flux_hi.col(face) - flux_lo.col(face)) / DX[face];
226 }
227 }
228
229 for (int p = 0; p < AMREX_SPACEDIM; ++p)
230 F(i, j, k, p) = f[p];
231 });
232 continue;
233 }
234
235 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
236
237 amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
238
239 Set::Vector f = Set::Vector::Zero();
240
241 Set::Vector u;
242 for (int p = 0; p < AMREX_SPACEDIM; p++) u(p) = U(i, j, k, p);
243
244
245 bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
246 AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
247
248 // Determine if a special stencil will be necessary for first derivatives
249 std::array<Numeric::StencilType, AMREX_SPACEDIM>
250 sten = Numeric::GetStencil(i, j, k, stencilbox);
251
252 // The displacement gradient tensor
253 Set::Matrix gradu; // gradu(i,j) = u_{i,j)
254
255 // Fill gradu
256 for (int p = 0; p < AMREX_SPACEDIM; p++)
257 {
258 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(U, i, j, k, p, DX, sten));,
259 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(U, i, j, k, p, DX, sten));,
260 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(U, i, j, k, p, DX, sten)););
261 }
262
263 Set::Scalar psi_avg = 1.0;
264 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
265
266 // Stress tensor computed using the model fab
267 Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
268
269 // Boundary conditions
270 /// \todo Important: we need a way to handle corners and edges.
271 amrex::IntVect m(AMREX_D_DECL(i, j, k));
272 if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
273 {
274 f = (*m_bc)(u, gradu, sig, i, j, k, stencilbox);
275 }
276 else
277 {
278
279
280 // The gradient of the displacement gradient tensor
281 // TODO - replace with this call. But not for this PR
282 //Set::Matrix3 gradgradu = Numeric::Hessian(U,i,j,k,DX,sten); // gradgradu[k](l,j) = u_{k,lj}
283 Set::Matrix3 gradgradu; // gradgradu[k](l,j) = u_{k,lj}
284
285 // Fill gradu and gradgradu
286 for (int p = 0; p < AMREX_SPACEDIM; p++)
287 {
288 // Diagonal terms:
289 AMREX_D_TERM(
290 gradgradu(p, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(U, i, j, k, p, DX));,
291 gradgradu(p, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(U, i, j, k, p, DX));,
292 gradgradu(p, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(U, i, j, k, p, DX)););
293
294 // Off-diagonal terms:
295 AMREX_D_TERM(
296 ,// 2D
297 gradgradu(p, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(U, i, j, k, p, DX));
298 gradgradu(p, 1, 0) = gradgradu(p, 0, 1);
299 ,// 3D
300 gradgradu(p, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(U, i, j, k, p, DX));
301 gradgradu(p, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(U, i, j, k, p, DX));
302 gradgradu(p, 2, 0) = gradgradu(p, 0, 2);
303 gradgradu(p, 2, 1) = gradgradu(p, 1, 2););
304 }
305
306 //
307 // Operator
308 //
309 // The return value is
310 // f = C(grad grad u) + grad(C)*grad(u)
311 // In index notation
312 // f_i = C_{ijkl,j} u_{k,l} + C_{ijkl}u_{k,lj}
313 //
314
315 f = (DDW(i, j, k) * gradgradu) * psi_avg;
316
317 if (!m_uniform)
318 {
319 MATRIX4
320 AMREX_D_DECL(Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(DDW, i, j, k, 0, DX, sten)),
321 Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(DDW, i, j, k, 0, DX, sten)),
322 Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(DDW, i, j, k, 0, DX, sten)));
323 f += (AMREX_D_TERM((Cgrad1 * gradu).col(0),
324 +(Cgrad2 * gradu).col(1),
325 +(Cgrad3 * gradu).col(2))) * (psi_avg);
326 }
327 if (m_psi_set)
328 {
329 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX);
330 gradpsi *= (1.0 - m_psi_small);
331 f += (DDW(i, j, k) * gradu) * gradpsi;
332 }
333
334 if (std::isnan(f(0)) || std::isnan(f(1)))
335 {
336 Util::Message(INFO," ================= ");
337 Util::Message(INFO,"amrlev=",amrlev);
338 Util::Message(INFO,"mglev=",mglev);
339 Util::Message(INFO,"i=",i," j=",j);
340 Util::Message(INFO,"f: ",f.transpose());
341 Util::Message(INFO,"U(i,j,k): ",U(i,j,k,0)," ",U(i,j,k,1));
342 Util::Message(INFO,"U(i-1,j,k): ",U(i-1,j,k,0)," ",U(i-1,j,k,1));
343 Util::Message(INFO,"U(i+1,j,k): ",U(i+1,j,k,0)," ",U(i-1,j,k,1));
344 Util::Message(INFO,"U(i,j-1,k): ",U(i,j-1,k,0)," ",U(i,j-1,k,1));
345 Util::Message(INFO,"U(i,j+1,k): ",U(i,j+1,k,0)," ",U(i,j+1,k,1));
346 Util::Message(INFO,"gradu: ",gradu);
347 Util::Message(INFO,"gradgradu[0]: ",gradgradu[0]);
348 Util::Message(INFO,"gradgradu[1]: ",gradgradu[1]);
349 Util::Message(INFO,"DDW (i ,j ): ",DDW(i,j,k));
350 Util::Message(INFO,"DDW (i-1,j ): ",DDW(i-1,j,k));
351 Util::Message(INFO,"DDW (i+1,j ): ",DDW(i+1,j,k));
352 Util::Message(INFO,"DDW (i ,j+1): ",DDW(i,j+1,k));
353 Util::Message(INFO,"DDW (i ,j-1): ",DDW(i,j-1,k));
354 Util::Message(INFO,"psi_av: ",psi_avg);
355 Util::Message(INFO,"psi_set: ",m_psi_set);
356 Util::Message(INFO," ================= ");
357 Util::Abort(INFO);
358 }
359
360 }
361 AMREX_D_TERM(F(i, j, k, 0) = f[0];, F(i, j, k, 1) = f[1];, F(i, j, k, 2) = f[2];);
362 });
363 }
364}
365
366
367
368template<int SYM>
369void
370Elastic<SYM>::Diagonal(int amrlev, int mglev, MultiFab& a_diag)
371{
372 BL_PROFILE("Operator::Elastic::Diagonal()");
373
374 // Conservative smoothing only consumes valid diagonal rows. Computing its
375 // ghost rows can cross into a neighboring coefficient FAB, where the local
376 // face data are not defined; FillBoundaryAndSync populates them below.
377 const amrex::IntVect diagonal_nghost = m_conservative_face_flux
378 ? amrex::IntVect::TheZeroVector() : a_diag.nGrowVect();
379 amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(
380 diagonal_nghost.max()));
381 domain.convert(amrex::IntVect::TheNodeVector());
382
383 amrex::Box stencilbox(m_geom[amrlev][mglev].growPeriodicDomain(
384 diagonal_nghost.max() + 1));
385 stencilbox.convert(amrex::IntVect::TheNodeVector());
386
387 const Real* DX = m_geom[amrlev][mglev].CellSize();
388
389 for (MFIter mfi(a_diag, false); mfi.isValid(); ++mfi)
390 {
391 Box bx = mfi.validbox().grow(diagonal_nghost) & domain;
392 amrex::Box tilebox = mfi.grownnodaltilebox() & bx;
393
394 amrex::Array4<MATRIX4> const& DDW = (*(m_ddw_mf[amrlev][mglev])).array(mfi);
395 amrex::Array4<Set::Scalar> const& diag = a_diag.array(mfi);
396 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][mglev]->array(mfi);
397
398 if (m_conservative_face_flux)
399 {
400 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
401 amrex::LoopConcurrentOnCpu(tilebox, [=] (int i, int j, int k)
402 {
403 const auto sten = Numeric::GetStencil(i, j, k, stencilbox);
404 const auto gradu =
405 Numeric::Gradient_Diagonal<Set::Matrix>(DX, sten);
406 Set::Scalar psi_avg = 1.0;
407 if (m_psi_set)
408 psi_avg = (1.0 - m_psi_small) *
409 Numeric::Interpolate::CellToNodeAverage(
410 psi, i, j, k, 0) + m_psi_small;
411 const int index[3] = {i, j, k};
412 const int lower[3] = {lo.x, lo.y, lo.z};
413 const int upper[3] = {hi.x, hi.y, hi.z};
414 bool on_boundary = false;
415 for (int dir = 0; dir < AMREX_SPACEDIM; ++dir)
416 on_boundary = on_boundary ||
417 index[dir] == lower[dir] || index[dir] == upper[dir];
418
419 for (int p = 0; p < AMREX_SPACEDIM; ++p)
420 {
421 diag(i, j, k, p) = 0.0;
422 if (on_boundary)
423 {
424 const Set::Matrix sig =
425 DDW(i, j, k) * gradu[p] * psi_avg;
426 Set::Vector u = Set::Vector::Zero();
427 u(p) = 1.0;
428 diag(i, j, k, p) =
429 (*m_bc)(u, gradu[p], sig, i, j, k, stencilbox)(p);
430 }
431 else
432 {
433 for (int face = 0; face < AMREX_SPACEDIM; ++face)
434 {
435 const int im = i - (face == 0);
436 const int jm = j - (face == 1);
437 const int km = k - (face == 2);
438 diag(i, j, k, p) -=
439 (DDW(i, j, k, face + 1)(p, face, p, face)
440 + DDW(im, jm, km, face + 1)(
441 p, face, p, face))
442 / (DX[face] * DX[face]);
443 }
444 }
445 }
446 });
447 continue;
448 }
449
450 const Dim3 lo = amrex::lbound(stencilbox), hi = amrex::ubound(stencilbox);
451
452 amrex::ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
453 {
454
455 // Determine if a special stencil will be necessary for first derivatives
456 std::array<Numeric::StencilType, AMREX_SPACEDIM>
457 sten = Numeric::GetStencil(i, j, k, stencilbox);
458
459 // gradu(i,j) = u_{i,j)
460 std::array<Set::Matrix,AMREX_SPACEDIM> gradu = Numeric::Gradient_Diagonal<Set::Matrix>(DX, sten);
461
462 // gradgradu[k](l,j) = u_{k,lj}
463 std::array<Set::Matrix3,AMREX_SPACEDIM> gradgradu = Numeric::Gradient_Diagonal<Set::Matrix3>(DX);
464
465
466 Set::Vector f = Set::Vector::Zero();
467
468 bool
469 AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z)),
470 AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
471
472 Set::Scalar psi_avg = 1.0;
473 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
474
475
476
477 for (int p = 0; p < AMREX_SPACEDIM; p++)
478 {
479
480 diag(i, j, k, p) = 0.0;
481
482
483 amrex::IntVect m(AMREX_D_DECL(i, j, k));
484 if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
485 {
486 Set::Matrix sig = DDW(i, j, k) * gradu[p] * psi_avg;
487 Set::Vector u = Set::Vector::Zero();
488 u(p) = 1.0;
489 f = (*m_bc)(u, gradu[p], sig, i, j, k, stencilbox);
490 diag(i, j, k, p) = f(p);
491 }
492 else
493 {
494 Set::Vector f = (DDW(i, j, k) * gradgradu[p]) * psi_avg;
495 diag(i, j, k, p) += f(p);
496 }
497
498#ifdef AMREX_DEBUG
499 if (std::isnan(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is nan at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
500 if (std::isinf(diag(i, j, k, p))) Util::Abort(INFO, "diagonal is inf at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
501 if (diag(i, j, k, p) == 0) Util::Abort(INFO, "diagonal is zero at (", i, ",", j, ",", k, "), amrlev=", amrlev, ", mglev=", mglev);
502#endif
503
504 }
505 });
506 }
507
508 a_diag.FillBoundaryAndSync(Geom(amrlev,mglev).periodicity());
509 nodalSync(amrlev,mglev,a_diag);
510}
511
512
513template<int SYM>
514void
515Elastic<SYM>::Error0x(int amrlev, int mglev, MultiFab& R0x, const MultiFab& x) const
516{
517 BL_PROFILE("Operator::Elastic::Error0x()");
518 Util::Message(INFO);
519
520 int ncomp = x.nComp();//getNComp();
521 int nghost = x.nGrow();
522
523 if (!m_diagonal_computed)
524 Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
525
526 amrex::MultiFab D0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
527 amrex::MultiFab AD0x(x.boxArray(), x.DistributionMap(), ncomp, nghost);
528
529 amrex::MultiFab::Copy(D0x, x, 0, 0, ncomp, nghost); // D0x = x
530 amrex::MultiFab::Divide(D0x, *m_diag[amrlev][mglev], 0, 0, ncomp, 0); // D0x = x/diag
531 amrex::MultiFab::Copy(AD0x, D0x, 0, 0, ncomp, nghost); // AD0x = D0x
532
533 Fapply(amrlev, mglev, AD0x, D0x); // AD0x = A * D0 * x
534
535 amrex::MultiFab::Copy(R0x, x, 0, 0, ncomp, nghost); // R0x = x
536 amrex::MultiFab::Subtract(R0x, AD0x, 0, 0, ncomp, nghost); // R0x = x - AD0x
537}
538
539
540template<int SYM>
541void
542Elastic<SYM>::FFlux(int /*amrlev*/, const MFIter& /*mfi*/,
543 const std::array<FArrayBox*, AMREX_SPACEDIM>& sigmafab,
544 const FArrayBox& /*ufab*/, const int /*face_only*/) const
545{
546 BL_PROFILE("Operator::Elastic::FFlux()");
547 Util::Message(INFO);
548 amrex::BaseFab<amrex::Real> AMREX_D_DECL(&fxfab = *sigmafab[0],
549 &fyfab = *sigmafab[1],
550 &fzfab = *sigmafab[2]);
551 AMREX_D_TERM(fxfab.setVal(0.0);,
552 fyfab.setVal(0.0);,
553 fzfab.setVal(0.0););
554
555}
556
557template<int SYM>
558void
559Elastic<SYM>::Strain(int amrlev,
560 amrex::MultiFab& a_eps,
561 const amrex::MultiFab& a_u,
562 bool voigt) const
563{
564 BL_PROFILE("Operator::Elastic::Strain()");
565
566 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
567 amrex::Box domain(m_geom[amrlev][0].Domain());
568 domain.convert(amrex::IntVect::TheNodeVector());
569
570
571 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
572 {
573 const Box& bx = mfi.tilebox();
574 amrex::Array4<amrex::Real> const& epsilon = a_eps.array(mfi);
575 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
576 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
577 {
578 Set::Matrix gradu;
579
580 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
581 = Numeric::GetStencil(i, j, k, domain);
582
583 // Fill gradu
584 for (int p = 0; p < AMREX_SPACEDIM; p++)
585 {
586 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
587 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
588 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
589 }
590
591 Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
592
593 if (voigt)
594 {
595 AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
596 ,
597 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(0, 1);
598 ,
599 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(1, 1); epsilon(i, j, k, 2) = eps(2, 2);
600 epsilon(i, j, k, 3) = eps(1, 2); epsilon(i, j, k, 4) = eps(2, 0); epsilon(i, j, k, 5) = eps(0, 1););
601 }
602 else
603 {
604 AMREX_D_PICK(epsilon(i, j, k, 0) = eps(0, 0);
605 ,
606 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1);
607 epsilon(i, j, k, 2) = eps(1, 0); epsilon(i, j, k, 3) = eps(1, 1);
608 ,
609 epsilon(i, j, k, 0) = eps(0, 0); epsilon(i, j, k, 1) = eps(0, 1); epsilon(i, j, k, 2) = eps(0, 2);
610 epsilon(i, j, k, 3) = eps(1, 0); epsilon(i, j, k, 4) = eps(1, 1); epsilon(i, j, k, 5) = eps(1, 2);
611 epsilon(i, j, k, 6) = eps(2, 0); epsilon(i, j, k, 7) = eps(2, 1); epsilon(i, j, k, 8) = eps(2, 2););
612 }
613 });
614 }
615}
616
617
618template<int SYM>
619void
620Elastic<SYM>::Stress(int amrlev,
621 amrex::MultiFab& a_sigma,
622 const amrex::MultiFab& a_u,
623 bool voigt, bool a_homogeneous)
624{
625 BL_PROFILE("Operator::Elastic::Stress()");
626 SetHomogeneous(a_homogeneous);
627
628 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
629 amrex::Box domain(m_geom[amrlev][0].Domain());
630 domain.convert(amrex::IntVect::TheNodeVector());
631
632 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
633 {
634 const Box& bx = mfi.tilebox();
635 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
636 amrex::Array4<amrex::Real> const& sigma = a_sigma.array(mfi);
637 amrex::Array4<Set::Scalar> const& psi = m_psi_mf[amrlev][0]->array(mfi);
638 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
639 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
640 {
641 Set::Matrix gradu;
642
643 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
644 = Numeric::GetStencil(i, j, k, domain);
645
646 // Fill gradu
647 for (int p = 0; p < AMREX_SPACEDIM; p++)
648 {
649 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
650 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
651 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
652 }
653
654 Set::Scalar psi_avg = 1.0;
655 if (m_psi_set) psi_avg = (1.0 - m_psi_small) * Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0) + m_psi_small;
656 Set::Matrix sig = (DDW(i, j, k) * gradu) * psi_avg;
657
658 if (voigt)
659 {
660 AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
661 ,
662 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(0, 1);
663 ,
664 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(1, 1); sigma(i, j, k, 2) = sig(2, 2);
665 sigma(i, j, k, 3) = sig(1, 2); sigma(i, j, k, 4) = sig(2, 0); sigma(i, j, k, 5) = sig(0, 1););
666 }
667 else
668 {
669 AMREX_D_PICK(sigma(i, j, k, 0) = sig(0, 0);
670 ,
671 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1);
672 sigma(i, j, k, 2) = sig(1, 0); sigma(i, j, k, 3) = sig(1, 1);
673 ,
674 sigma(i, j, k, 0) = sig(0, 0); sigma(i, j, k, 1) = sig(0, 1); sigma(i, j, k, 2) = sig(0, 2);
675 sigma(i, j, k, 3) = sig(1, 0); sigma(i, j, k, 4) = sig(1, 1); sigma(i, j, k, 5) = sig(1, 2);
676 sigma(i, j, k, 6) = sig(2, 0); sigma(i, j, k, 7) = sig(2, 1); sigma(i, j, k, 8) = sig(2, 2););
677 }
678 });
679 }
680}
681
682
683template<int SYM>
684void
685Elastic<SYM>::Energy(int amrlev,
686 amrex::MultiFab& a_energy,
687 const amrex::MultiFab& a_u, bool a_homogeneous)
688{
689 BL_PROFILE("Operator::Elastic::Energy()");
690 SetHomogeneous(a_homogeneous);
691
692 amrex::Box domain(m_geom[amrlev][0].Domain());
693 domain.convert(amrex::IntVect::TheNodeVector());
694
695 const amrex::Real* DX = m_geom[amrlev][0].CellSize();
696
697 for (MFIter mfi(a_u, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
698 {
699 const Box& bx = mfi.tilebox();
700 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& DDW = (*(m_ddw_mf[amrlev][0])).array(mfi);
701 amrex::Array4<amrex::Real> const& energy = a_energy.array(mfi);
702 amrex::Array4<const amrex::Real> const& u = a_u.array(mfi);
703 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
704 {
705 Set::Matrix gradu;
706
707 std::array<Numeric::StencilType, AMREX_SPACEDIM> sten
708 = Numeric::GetStencil(i, j, k, domain);
709
710 // Fill gradu
711 for (int p = 0; p < AMREX_SPACEDIM; p++)
712 {
713 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
714 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
715 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
716 }
717
718 Set::Matrix eps = .5 * (gradu + gradu.transpose());
719 Set::Matrix sig = DDW(i, j, k) * gradu;
720
721 // energy(i,j,k) = (gradu.transpose() * sig).trace();
722
723 //Util::Abort(INFO,"Fix this"); //
724 //energy(i,j,k) = C(i,j,k).W(gradu);
725 for (int m = 0; m < AMREX_SPACEDIM; m++)
726 {
727 for (int n = 0; n < AMREX_SPACEDIM; n++)
728 {
729 energy(i, j, k) += .5 * sig(m, n) * eps(m, n);
730 }
731 }
732 });
733 }
734}
735
736template<int SYM>
737void
738Elastic<SYM>::averageDownCoeffs()
739{
740 BL_PROFILE("Elastic::averageDownCoeffs()");
741
742 if (m_average_down_coeffs)
743 for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
744 averageDownCoeffsDifferentAmrLevels(amrlev);
745
746 averageDownCoeffsSameAmrLevel(0);
747 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
748 {
749 for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
750 {
751 if (m_ddw_mf[amrlev][mglev]) {
752 FillBoundaryCoeff(*m_ddw_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
753 FillBoundaryCoeff(*m_psi_mf[amrlev][mglev], Geom(amrlev,mglev).periodicity());
754 }
755 }
756 }
757}
758
759template<int SYM>
760void
761Elastic<SYM>::averageDownCoeffsDifferentAmrLevels(int fine_amrlev)
762{
763 BL_PROFILE("Operator::Elastic::averageDownCoeffsDifferentAmrLevels()");
764 Util::Assert(INFO, TEST(fine_amrlev > 0));
765
766 const int crse_amrlev = fine_amrlev - 1;
767
768 MultiTab& crse_ddw = *m_ddw_mf[crse_amrlev][0];
769 MultiTab& fine_ddw = *m_ddw_mf[fine_amrlev][0];
770 const int ncomp = crse_ddw.nComp();
771
772 amrex::Box cdomain(m_geom[crse_amrlev][0].Domain());
773 cdomain.convert(amrex::IntVect::TheNodeVector());
774
775 const Geometry& cgeom = m_geom[crse_amrlev][0];
776
777 const BoxArray& fba = fine_ddw.boxArray();
778 const DistributionMapping& fdm = fine_ddw.DistributionMap();
779
780 MultiTab fine_ddw_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
781 fine_ddw_for_coarse.ParallelCopy(crse_ddw, 0, 0, ncomp, 0, 0, cgeom.periodicity());
782
783 const int coarse_fine_node = 1;
784 const int fine_fine_node = 2;
785
786 amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
787 nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
788
789 amrex::iMultiFab cellmask(amrex::convert(amrex::coarsen(fba, 2), amrex::IntVect::TheCellVector()), fdm, 1, 2);
790 cellmask.ParallelCopy(*m_cc_fine_mask[crse_amrlev], 0, 0, 1, 1, 1, cgeom.periodicity());
791
792 for (MFIter mfi(fine_ddw_for_coarse, false); mfi.isValid(); ++mfi)
793 {
794 const Box& bx = mfi.validbox();
795
796 amrex::Array4<const int> const& nmask = nodemask.array(mfi);
797 //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
798
799 amrex::Array4<MATRIX4> const& cdata = fine_ddw_for_coarse.array(mfi);
800 amrex::Array4<const MATRIX4> const& fdata = fine_ddw.array(mfi);
801
802 const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
803
804 for (int n = 0; n < ncomp; n++)
805 {
806 // I,J,K == coarse coordinates
807 // i,j,k == fine coordinates
808 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
809 int i = I * 2, j = J * 2, k = K * 2;
810
811 if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
812 {
813 if (n > 0)
814 {
815 const int face = n - 1;
816 cdata(I, J, K, n) = 0.5 * (
817 fdata(i, j, k, n)
818 + fdata(i + (face == 0), j + (face == 1),
819 k + (face == 2), n));
820 return;
821 }
822 if ((I == lo.x || I == hi.x) &&
823 (J == lo.y || J == hi.y) &&
824 (K == lo.z || K == hi.z)) // Corner
825 cdata(I, J, K, n) = fdata(i, j, k, n);
826 else if ((J == lo.y || J == hi.y) &&
827 (K == lo.z || K == hi.z)) // X edge
828 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;
829 else if ((K == lo.z || K == hi.z) &&
830 (I == lo.x || I == hi.x)) // Y edge
831 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;
832 else if ((I == lo.x || I == hi.x) &&
833 (J == lo.y || J == hi.y)) // Z edge
834 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;
835 else if (I == lo.x || I == hi.x) // X face
836 cdata(I, J, K, n) =
837 (fdata(i, j - 1, k - 1, n) + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j + 1, k - 1, n)
838 + fdata(i, j - 1, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j + 1, k, n) * 2.0
839 + fdata(i, j - 1, k + 1, n) + fdata(i, j, k + 1, n) * 2.0 + fdata(i, j + 1, k + 1, n)) / 16.0;
840 else if (J == lo.y || J == hi.y) // Y face
841 cdata(I, J, K, n) =
842 (fdata(i - 1, j, k - 1, n) + fdata(i - 1, j, k, n) * 2.0 + fdata(i - 1, j, k + 1, n)
843 + fdata(i, j, k - 1, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i, j, k + 1, n) * 2.0
844 + fdata(i + 1, j, k - 1, n) + fdata(i + 1, j, k, n) * 2.0 + fdata(i + 1, j, k + 1, n)) / 16.0;
845 else if (K == lo.z || K == hi.z) // Z face
846 cdata(I, J, K, n) =
847 (fdata(i - 1, j - 1, k, n) + fdata(i, j - 1, k, n) * 2.0 + fdata(i + 1, j - 1, k, n)
848 + fdata(i - 1, j, k, n) * 2.0 + fdata(i, j, k, n) * 4.0 + fdata(i + 1, j, k, n) * 2.0
849 + fdata(i - 1, j + 1, k, n) + fdata(i, j + 1, k, n) * 2.0 + fdata(i + 1, j + 1, k, n)) / 16.0;
850 else // Interior
851 cdata(I, J, K, n) =
852 (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) +
853 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
854 +
855 (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) +
856 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) +
857 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
858 +
859 (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
860 fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
861 +
862 fdata(i, j, k, n) / 8.0;
863
864#ifdef AMREX_DEBUG
865 if (cdata(I, J, K, n).contains_nan()) Util::Abort(INFO, "restricted model is nan at (", i, ",", j, ",", k, "), fine_amrlev=", fine_amrlev);
866#endif
867 }
868
869 });
870 }
871 }
872
873 // Copy the fine residual restricted onto the coarse grid
874 // into the final residual.
875
876 crse_ddw.ParallelCopy(fine_ddw_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
877 //const int mglev = 0;
878 //Util::RealFillBoundary(crse_ddw, m_geom[crse_amrlev][mglev]);
879
880 FillBoundaryCoeff(crse_ddw,Geom(fine_amrlev,0).periodicity());
881
882}
883
884
885
886template<int SYM>
887void
888Elastic<SYM>::averageDownCoeffsSameAmrLevel(int amrlev)
889{
890 BL_PROFILE("Elastic::averageDownCoeffsSameAmrLevel()");
891
892 for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev)
893 {
894 amrex::Box cdomain(m_geom[amrlev][mglev].growPeriodicDomain(2));
895 cdomain.convert(amrex::IntVect::TheNodeVector());
896 amrex::Box fdomain(m_geom[amrlev][mglev - 1].Domain());
897 fdomain.convert(amrex::IntVect::TheNodeVector());
898
899 MultiTab& crse = *m_ddw_mf[amrlev][mglev];
900 MultiTab& fine = *m_ddw_mf[amrlev][mglev - 1];
901 const int ncomp = crse.nComp();
902
903 amrex::BoxArray crseba = crse.boxArray();
904 amrex::BoxArray fineba = fine.boxArray();
905
906 BoxArray newba = crseba;
907 newba.refine(2);
908 MultiTab fine_on_crseba;
909 fine_on_crseba.define(newba, crse.DistributionMap(), ncomp, 4);
910 fine_on_crseba.ParallelCopy(fine, 0, 0, ncomp, 2, 4,
911 m_geom[amrlev][mglev-1].periodicity());
912 /* ine_on_crseba.FillBoundaryAndSync(m_geom[amrlev][mglev-1].periodicity()); */
913
914 for (MFIter mfi(crse, false); mfi.isValid(); ++mfi)
915 {
916
917 Box bx = mfi.grownnodaltilebox() & cdomain;
918 /*Box bx = mfi.grownnodaltilebox(-1,1) & cdomain;*/
919
920 amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>> const& fdata = fine_on_crseba.array(mfi);
921 amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, SYM>> const& cdata = crse.array(mfi);
922
923 const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
924 /*const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);*/
925
926 // I,J,K == coarse coordinates
927 // i,j,k == fine coordinates
928 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
929 int i = 2 * I, j = 2 * J, k = 2 * K;
930
931 if ((I == lo.x || I == hi.x) &&
932 (J == lo.y || J == hi.y) &&
933 (K == lo.z || K == hi.z)) // Corner
934 cdata(I, J, K) = fdata(i, j, k);
935 else if ((J == lo.y || J == hi.y) &&
936 (K == lo.z || K == hi.z)) // X edge
937 cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
938 else if ((K == lo.z || K == hi.z) &&
939 (I == lo.x || I == hi.x)) // Y edge
940 cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
941 else if ((I == lo.x || I == hi.x) &&
942 (J == lo.y || J == hi.y)) // Z edge
943 cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
944 else if (I == lo.x || I == hi.x) // X face
945 cdata(I, J, K) =
946 (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
947 + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
948 + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
949 else if (J == lo.y || J == hi.y) // Y face
950 cdata(I, J, K) =
951 (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
952 + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
953 + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
954 else if (K == lo.z || K == hi.z) // Z face
955 cdata(I, J, K) =
956 (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
957 + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
958 + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
959 else // Interior
960 cdata(I, J, K) =
961 (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) +
962 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
963 +
964 (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
965 fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
966 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
967 +
968 (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
969 fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
970 +
971 fdata(i, j, k) / 8.0;
972
973#ifdef AMREX_DEBUG
974 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);
975#endif
976 });
977
978 for (int n = 1; n < ncomp; ++n)
979 {
980 const int face = n - 1;
981 amrex::LoopConcurrentOnCpu(bx, [=] (int I, int J, int K)
982 {
983 const int i = 2 * I, j = 2 * J, k = 2 * K;
984 cdata(I, J, K, n) = 0.5 * (
985 fdata(i, j, k, n)
986 + fdata(i + (face == 0), j + (face == 1),
987 k + (face == 2), n));
988 });
989 }
990 }
991 FillBoundaryCoeff(crse, Geom(amrlev,mglev).periodicity());
992
993
994 if (!m_psi_set) continue;
995
996 amrex::Box cdomain_cell(m_geom[amrlev][mglev].Domain());
997 amrex::Box fdomain_cell(m_geom[amrlev][mglev - 1].Domain());
998 MultiFab& crse_psi = *m_psi_mf[amrlev][mglev];
999 MultiFab& fine_psi = *m_psi_mf[amrlev][mglev - 1];
1000 MultiFab fine_psi_on_crseba;
1001 fine_psi_on_crseba.define(newba.convert(amrex::IntVect::TheCellVector()), crse_psi.DistributionMap(), 1, 1);
1002 fine_psi_on_crseba.ParallelCopy(fine_psi, 0, 0, 1, 1, 1, m_geom[amrlev][mglev].periodicity());
1003
1004 for (MFIter mfi(crse_psi, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
1005 {
1006 Box bx = mfi.tilebox();
1007 bx = bx & cdomain_cell;
1008
1009 amrex::Array4<const Set::Scalar> const& fdata = fine_psi_on_crseba.array(mfi);
1010 amrex::Array4<Set::Scalar> const& cdata = crse_psi.array(mfi);
1011
1012 const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
1013
1014 // I,J,K == coarse coordinates
1015 // i,j,k == fine coordinates
1016 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
1017 int i = 2 * I, j = 2 * J, k = 2 * K;
1018
1019 if ((I == lo.x || I == hi.x) &&
1020 (J == lo.y || J == hi.y) &&
1021 (K == lo.z || K == hi.z)) // Corner
1022 cdata(I, J, K) = fdata(i, j, k);
1023 else if ((J == lo.y || J == hi.y) &&
1024 (K == lo.z || K == hi.z)) // X edge
1025 cdata(I, J, K) = fdata(i - 1, j, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i + 1, j, k) * 0.25;
1026 else if ((K == lo.z || K == hi.z) &&
1027 (I == lo.x || I == hi.x)) // Y edge
1028 cdata(I, J, K) = fdata(i, j - 1, k) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j + 1, k) * 0.25;
1029 else if ((I == lo.x || I == hi.x) &&
1030 (J == lo.y || J == hi.y)) // Z edge
1031 cdata(I, J, K) = fdata(i, j, k - 1) * 0.25 + fdata(i, j, k) * 0.5 + fdata(i, j, k + 1) * 0.25;
1032 else if (I == lo.x || I == hi.x) // X face
1033 cdata(I, J, K) =
1034 (fdata(i, j - 1, k - 1) + fdata(i, j, k - 1) * 2.0 + fdata(i, j + 1, k - 1)
1035 + fdata(i, j - 1, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j + 1, k) * 2.0
1036 + fdata(i, j - 1, k + 1) + fdata(i, j, k + 1) * 2.0 + fdata(i, j + 1, k + 1)) / 16.0;
1037 else if (J == lo.y || J == hi.y) // Y face
1038 cdata(I, J, K) =
1039 (fdata(i - 1, j, k - 1) + fdata(i - 1, j, k) * 2.0 + fdata(i - 1, j, k + 1)
1040 + fdata(i, j, k - 1) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i, j, k + 1) * 2.0
1041 + fdata(i + 1, j, k - 1) + fdata(i + 1, j, k) * 2.0 + fdata(i + 1, j, k + 1)) / 16.0;
1042 else if (K == lo.z || K == hi.z) // Z face
1043 cdata(I, J, K) =
1044 (fdata(i - 1, j - 1, k) + fdata(i, j - 1, k) * 2.0 + fdata(i + 1, j - 1, k)
1045 + fdata(i - 1, j, k) * 2.0 + fdata(i, j, k) * 4.0 + fdata(i + 1, j, k) * 2.0
1046 + fdata(i - 1, j + 1, k) + fdata(i, j + 1, k) * 2.0 + fdata(i + 1, j + 1, k)) / 16.0;
1047 else // Interior
1048 cdata(I, J, K) =
1049 (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) +
1050 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
1051 +
1052 (fdata(i, j - 1, k - 1) + fdata(i, j - 1, k + 1) + fdata(i, j + 1, k - 1) + fdata(i, j + 1, k + 1) +
1053 fdata(i - 1, j, k - 1) + fdata(i + 1, j, k - 1) + fdata(i - 1, j, k + 1) + fdata(i + 1, j, k + 1) +
1054 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
1055 +
1056 (fdata(i - 1, j, k) + fdata(i, j - 1, k) + fdata(i, j, k - 1) +
1057 fdata(i + 1, j, k) + fdata(i, j + 1, k) + fdata(i, j, k + 1)) / 16.0
1058 +
1059 fdata(i, j, k) / 8.0;
1060 });
1061 }
1062 FillBoundaryCoeff(crse_psi, Geom(amrlev,mglev).periodicity());
1063
1064 }
1065}
1066
1067template<int SYM>
1068void
1069Elastic<SYM>::FillBoundaryCoeff(MultiTab& sigma, const amrex::Periodicity& p)
1070{
1071 sigma.setMultiGhost(true);
1072 sigma.FillBoundaryAndSync(p);
1073}
1074
1075template<int SYM>
1076void
1077Elastic<SYM>::FillBoundaryCoeff(MultiFab& psi, const amrex::Periodicity& p)
1078{
1079 psi.setMultiGhost(true);
1080 psi.FillBoundaryAndSync(p);
1081}
1082
1083template class Elastic<Set::Sym::Major>;
1084template class Elastic<Set::Sym::Isotropic>;
1085template class Elastic<Set::Sym::MajorMinor>;
1086template class Elastic<Set::Sym::Diagonal>;
1087
1088}
TEST
#define TEST(x)
Definition Util.H:25
INFO
#define INFO
Definition Util.H:24
Operator::Elastic::SetModel
void SetModel(Set::Matrix4< AMREX_SPACEDIM, SYM > &a_model)
Definition Elastic.cpp:67
Operator::Elastic< T::sym >::MultiTab
amrex::FabArray< TArrayBox > MultiTab
Definition Elastic.H:26
Operator::Elastic::Stress
void Stress(int amrlev, amrex::MultiFab &sigmafab, const amrex::MultiFab &ufab, bool voigt=false, bool a_homogeneous=false)
Compute stress given the displacement field by.
Definition Elastic.cpp:620
Operator::Elastic::FFlux
virtual void FFlux(int amrlev, const MFIter &mfi, const std::array< FArrayBox *, AMREX_SPACEDIM > &flux, const FArrayBox &sol, const int face_only=0) const final
Definition Elastic.cpp:542
Operator::Elastic::Error0x
void Error0x(int amrlev, int mglev, MultiFab &R0x, const MultiFab &x) const
Definition Elastic.cpp:515
Operator::Elastic::define
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FArrayBox > const * > &a_factory={}, bool a_conservative_face_flux=false)
Definition Elastic.cpp:28
Operator::Elastic::FillBoundaryCoeff
void FillBoundaryCoeff(MultiTab &sigma, const amrex::Periodicity &p)
Definition Elastic.cpp:1069
Operator::Elastic::averageDownCoeffsSameAmrLevel
void averageDownCoeffsSameAmrLevel(int amrlev)
Update coarse-level AMR coefficients with data from fine level.
Definition Elastic.cpp:888
Operator::Elastic::~Elastic
virtual ~Elastic()
Definition Elastic.cpp:23
Operator::Elastic::Strain
void Strain(int amrlev, amrex::MultiFab &epsfab, const amrex::MultiFab &ufab, bool voigt=false) const
Compute strain given the displacement field by.
Definition Elastic.cpp:559
Operator::Elastic::Energy
void Energy(int amrlev, amrex::MultiFab &energy, const amrex::MultiFab &u, bool a_homogeneous=false)
Compute energy density given the displacement field by.
Definition Elastic.cpp:685
Operator::Elastic::averageDownCoeffsDifferentAmrLevels
void averageDownCoeffsDifferentAmrLevels(int fine_amrlev)
Definition Elastic.cpp:761
Operator::Elastic::Fapply
virtual void Fapply(int amrlev, int mglev, MultiFab &out, const MultiFab &in) const override final
Definition Elastic.cpp:158
Operator::Elastic::Diagonal
virtual void Diagonal(int amrlev, int mglev, amrex::MultiFab &diag) override
Definition Elastic.cpp:370
Operator::Elastic::averageDownCoeffs
virtual void averageDownCoeffs() override
Definition Elastic.cpp:738
Operator::Elastic::SetPsi
void SetPsi(int amrlev, const amrex::MultiFab &a_psi)
Definition Elastic.cpp:137
Numeric::CellGradientOnNode
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Vector CellGradientOnNode(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:705
Numeric::FaceGradient
Set::Matrix FaceGradient(const amrex::Array4< const Set::Vector > &f, const int i, const int j, const int k, const int face, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:796
Numeric::GetStencil
static AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition Stencil.H:51
Numeric::Gradient
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType())
Definition Stencil.H:687
Operator
Documentation for operator namespace.
Definition Diagonal.cpp:14
Operator::AMREX_D_DECL
constexpr amrex::IntVect AMREX_D_DECL(Operator< Grid::Cell >::dx, Operator< Grid::Cell >::dy, Operator< Grid::Cell >::dz)
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:21
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:24
Util::Assert
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE void Assert(const char *file, const char *func, int line, const char *smt, bool pass, Args const &... args)
Definition Util.H:60
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:140
Util::Abort
AMREX_GPU_HOST_DEVICE void Abort(const char *msg)
Definition Util.cpp:406
Numeric::Interpolate::CellToNodeAverage
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType())
Definition Stencil.H:1445
Generated by doxygen 1.9.8