OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
AmrTrilinearInterpolater.hpp
Go to the documentation of this file.
1 //
2 // Class AmrTrilinearInterpolater
3 // Trilinear interpolation of data on coarse cells to fine cells.
4 //
5 // Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
6 // All rights reserved
7 //
8 // Implemented as part of the PhD thesis
9 // "Precise Simulations of Multibunches in High Intensity Cyclotrons"
10 //
11 // This file is part of OPAL.
12 //
13 // OPAL is free software: you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation, either version 3 of the License, or
16 // (at your option) any later version.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
20 //
21 template <class Level>
23  : AmrInterpolater<Level>(2 << (AMREX_SPACEDIM - 1))
24 { }
25 
26 
27 template <class Level>
29  const AmrIntVect_t& iv,
30  const basefab_t& fab,
31  umap_t& map,
32  const scalar_t& scale,
33  Level* mglevel)
34 {
35  /* lower left coarse cell (i, j, k)
36  * floor( i - 0.5 ) / rr[0]
37  * floor( j - 0.5 ) / rr[1]
38  * floor( k - 0.5 ) / rr[2]
39  */
40  AmrIntVect_t civ;
41  for (int d = 0; d < AMREX_SPACEDIM; ++d) {
42 
43  scalar_t tmp = iv[d] - 0.5;
44  if ( std::signbit(tmp) )
45  civ[d] = std::floor(tmp);
46  else
47  civ[d] = tmp;
48  }
49 
50  civ.coarsen(mglevel->refinement());
51 
52  // ref ratio 2 only
53  scalar_t dx = 0.5 * ( iv[0] - civ[0] * 2 ) - 0.25;
54  scalar_t dy = 0.5 * ( iv[1] - civ[1] * 2 ) - 0.25;
55 #if AMREX_SPACEDIM == 3
56  scalar_t dz = 0.5 * ( iv[2] - civ[2] * 2 ) - 0.25;
57 #endif
58 
59  scalar_t xdiff = 1.0 - dx;
60  scalar_t ydiff = 1.0 - dy;
61 #if AMREX_SPACEDIM == 3
62  scalar_t zdiff = 1.0 - dz;
63 #endif
64  // (i, j, k)
65  go_t crse_gidx = mglevel->serialize(civ);
66  scalar_t value = AMREX_D_TERM(xdiff, * ydiff, * zdiff) * scale;
67 
68  if ( !mglevel->applyBoundary(civ, fab, map, value) )
69  map[crse_gidx] += value;
70 
71  // (i+1, j, k)
72  AmrIntVect_t tmp(D_DECL(civ[0]+1, civ[1], civ[2]));
73  value = AMREX_D_TERM(dx, * ydiff, * zdiff) * scale;
74  if ( !mglevel->applyBoundary(tmp, map, value) )
75  map[mglevel->serialize(tmp)] += value;
76 
77  // (i, j+1, k)
78  tmp = AmrIntVect_t(D_DECL(civ[0], civ[1]+1, civ[2]));
79  value = AMREX_D_TERM(xdiff, * dy, * zdiff) * scale;
80  if ( !mglevel->applyBoundary(tmp, map, value) )
81  map[mglevel->serialize(tmp)] += value;
82 
83  // (i+1, j+1, k)
84  tmp = AmrIntVect_t(D_DECL(civ[0]+1, civ[1]+1, civ[2]));
85  value = AMREX_D_TERM(dx, * dy, * zdiff) * scale;
86  if ( !mglevel->applyBoundary(tmp, map, value) )
87  map[mglevel->serialize(tmp)] += value;
88 
89 #if AMREX_SPACEDIM == 3
90  // (i, j, k+1)
91  tmp = AmrIntVect_t(D_DECL(civ[0], civ[1], civ[2]+1));
92  value = AMREX_D_TERM(xdiff, * ydiff, * dz) * scale;
93  if ( !mglevel->applyBoundary(tmp, map, value) )
94  map[mglevel->serialize(tmp)] += value;
95 
96  // (i+1, j, k+1)
97  tmp = AmrIntVect_t(D_DECL(civ[0]+1, civ[1], civ[2]+1));
98  value = AMREX_D_TERM(dx, * ydiff, * dz) * scale;
99  if ( !mglevel->applyBoundary(tmp, map, value) )
100  map[mglevel->serialize(tmp)] += value;
101 
102  // (i, j+1, k+1)
103  tmp = AmrIntVect_t(D_DECL(civ[0], civ[1]+1, civ[2]+1));
104  value = AMREX_D_TERM(xdiff, * dy, * dz) * scale;
105  if ( !mglevel->applyBoundary(tmp, map, value) )
106  map[mglevel->serialize(tmp)] += value;
107 
108  // (i+1, j+1, k+1)
109  tmp = AmrIntVect_t(D_DECL(civ[0]+1, civ[1]+1, civ[2]+1));
110  value = AMREX_D_TERM(dx, * dy, * dz) * scale;
111  if ( !mglevel->applyBoundary(tmp, map, value) )
112  map[mglevel->serialize(tmp)] += value;
113 #endif
114 }
115 
116 
117 template <class Level>
119  const AmrIntVect_t& /*iv*/,
120  umap_t& /*map*/,
121  const scalar_t& /*scale*/,
122  lo_t /*dir*/, lo_t /*shift*/, const basefab_t& /*rfab*/,
123  const AmrIntVect_t& /*riv*/,
124  Level* /*mglevel*/)
125 {
126  // do nothing
127 }
128 
129 
130 template <class Level>
132  const AmrIntVect_t& iv,
133  umap_t& map,
134  const scalar_t& scale,
135  lo_t /*dir*/, lo_t /*shift*/, const basefab_t& fab,
136  Level* mglevel)
137 {
138  /*
139  * The AmrTrilinearInterpolater interpolates directly to the
140  * fine ghost cell.
141  */
142  this->stencil(iv, fab, map, scale, mglevel);
143 }
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:733
amrex::IntVect AmrIntVect_t
Definition: AmrDefs.h:48
amr::AmrIntVect_t AmrIntVect_t
Level::scalar_t scalar_t
Level::lo_t lo_t
Level::umap_t umap_t
Level::go_t go_t
Level::basefab_t basefab_t
void stencil(const AmrIntVect_t &iv, const basefab_t &fab, umap_t &map, const scalar_t &scale, Level *mglevel)
void fine(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, const basefab_t &fab, Level *mglevel)
void coarse(const AmrIntVect_t &iv, umap_t &map, const scalar_t &scale, lo_t dir, lo_t shift, const basefab_t &rfab, const AmrIntVect_t &riv, Level *mglevel)