AmrMultiGridLevel.hpp 8.33 KB
Newer Older
1 2 3 4
//
// Class AmrMultiGridLevel
//   This class represents a single AMR level, i.e. it stores all matrices
//   and vectors of a level.
5 6 7 8
//
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
9 10 11 12 13 14 15 16 17 18 19 20 21
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
frey_m's avatar
frey_m committed
22 23 24
#define AMR_NO_SCALE false


25 26
template <class MatrixType, class VectorType>
AmrMultiGridLevel<MatrixType,
frey_m's avatar
frey_m committed
27 28
                  VectorType>::AmrMultiGridLevel(const Vector_t& meshScaling,
                                                 const amrex::BoxArray& _grids,
29 30 31
                                                 const amrex::DistributionMapping& _dmap,
                                                 const AmrGeometry_t& _geom,
                                                 const AmrIntVect_t& rr,
32
                                                 const boundary_t* bc,
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
                                                 const Teuchos::RCP<comm_t>& comm,
                                                 const Teuchos::RCP<node_t>& node)
    : grids(_grids),
      dmap(_dmap),
      geom(_geom),
      map_p(Teuchos::null),
      Anf_p(Teuchos::null),
      R_p(Teuchos::null),
      I_p(Teuchos::null),
      Bcrse_p(Teuchos::null),
      Bfine_p(Teuchos::null),
      Awf_p(Teuchos::null),
      rho_p(Teuchos::null),
      phi_p(Teuchos::null),
      residual_p(Teuchos::null),
      error_p(Teuchos::null),
      UnCovered_p(Teuchos::null),
50
      refmask(nullptr),
frey_m's avatar
frey_m committed
51 52
      crsemask(nullptr),
      rr_m(rr)
53 54 55 56 57
{
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        G_p[j] = Teuchos::null;
        
        nr_m[j] = _geom.Domain().length(j);
58
        
frey_m's avatar
frey_m committed
59 60 61 62 63 64 65 66 67 68
#if AMR_NO_SCALE
        // mesh spacing in particle rest frame
        dx_m[j] = geom.CellSize(j);
        invdx_m[j] = geom.InvCellSize(j);
#else
        // mesh spacing in particle rest frame
        dx_m[j] = meshScaling[j] * geom.CellSize(j);
        invdx_m[j] = meshScaling[j] * geom.InvCellSize(j);
#endif
        
69
        bc_mp[j] = bc[j];
70 71
    }
    
frey_m's avatar
frey_m committed
72
    this->buildLevelMask();
73
    
frey_m's avatar
frey_m committed
74
    this->buildMap(comm, node);
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
    
    
    residual_p = Teuchos::rcp( new vector_t(map_p, false) );
    error_p = Teuchos::rcp( new vector_t(map_p, false) );
}


template <class MatrixType, class VectorType>
AmrMultiGridLevel<MatrixType, VectorType>::~AmrMultiGridLevel()
{
    map_p = Teuchos::null;
    
    Anf_p = Teuchos::null;
    R_p = Teuchos::null;
    I_p = Teuchos::null;
    Bcrse_p = Teuchos::null;
    Bfine_p = Teuchos::null;
    Awf_p = Teuchos::null;
    
    for (int j = 0; j < AMREX_SPACEDIM; ++j)
        G_p[j] = Teuchos::null;
    
    UnCovered_p = Teuchos::null;
    
    rho_p = Teuchos::null;
    phi_p = Teuchos::null;
    residual_p = Teuchos::null;
    error_p = Teuchos::null;
}


template <class MatrixType, class VectorType>
frey_m's avatar
frey_m committed
107 108
typename AmrMultiGridLevel<MatrixType, VectorType>::go_t
AmrMultiGridLevel<MatrixType, VectorType>::serialize(const AmrIntVect_t& iv) const {
109
#if AMREX_SPACEDIM == 3
110
    return iv[0] + (iv[1] + nr_m[1] * iv[2]) * nr_m[0];
111
#else
112
    return iv[0] + iv[1] * nr_m[0];
113 114 115 116 117 118
#endif
}


template <class MatrixType, class VectorType>
bool AmrMultiGridLevel<MatrixType, VectorType>::isBoundary(const AmrIntVect_t& iv) const {
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
    // it doesn't matter with which direction we check, since it checks all
    return bc_mp[0]->isBoundary(iv, &nr_m[0]);
}


template <class MatrixType, class VectorType>
bool AmrMultiGridLevel<MatrixType, VectorType>::applyBoundary(const AmrIntVect_t& iv,
                                                              umap_t& map,
                                                              const scalar_t& value)
{
    bool applied = false;
    for (int d = 0; d < AMREX_SPACEDIM; ++d) {
        if ( bc_mp[d]->isBoundary(iv, d, &nr_m[0]) ) {
            applied = true;
            bc_mp[d]->apply(iv, d, map, value, this, &nr_m[0]);
        }
    }
    return applied;
}


template <class MatrixType, class VectorType>
bool AmrMultiGridLevel<MatrixType, VectorType>::applyBoundary(const AmrIntVect_t& iv,
                                                              const basefab_t& fab,
                                                              umap_t& map,
                                                              const scalar_t& value)
{
    if ( fab(iv) != Mask::PHYSBNDRY )
        return false;
    
    bool applied = false;
    for (int d = 0; d < AMREX_SPACEDIM; ++d) {
        if ( bc_mp[d]->isBoundary(iv, d, &nr_m[0]) ) {
            applied = true;
            bc_mp[d]->apply(iv, d, map, value, this, &nr_m[0]);
        }
    }
    return applied;
157 158 159 160 161
}


template <class MatrixType, class VectorType>
void AmrMultiGridLevel<MatrixType, VectorType>::applyBoundary(const AmrIntVect_t& iv,
162 163 164
                                                              const lo_t& dir,
                                                              umap_t& map,
                                                              const scalar_t& value)
165
{
166
    bc_mp[dir]->apply(iv, dir, map, value, this, &nr_m[0]);
167 168 169 170
}


template <class MatrixType, class VectorType>
frey_m's avatar
frey_m committed
171
void AmrMultiGridLevel<MatrixType, VectorType>::buildLevelMask() {
172
    amrex::Periodicity period(AmrIntVect_t(D_DECL(0, 0, 0)));
173
    mask.reset(new mask_t(grids, dmap, 1, 1));
174
    mask->BuildMask(geom.Domain(), period,
175 176
                    Mask::COVERED, Mask::BNDRY,
                    Mask::PHYSBNDRY, Mask::INTERIOR);
frey_m's avatar
frey_m committed
177
    mask->FillBoundary(period);
178 179 180 181
}


template <class MatrixType, class VectorType>
182
const amr::AmrIntVect_t& AmrMultiGridLevel<MatrixType, VectorType>::refinement() const {
183 184 185 186
    return rr_m;
}


frey_m's avatar
frey_m committed
187 188 189 190 191 192 193
template <class MatrixType, class VectorType>
const amr::scalar_t* AmrMultiGridLevel<MatrixType, VectorType>::cellSize() const {
    return dx_m;
}


template <class MatrixType, class VectorType>
frey_m's avatar
frey_m committed
194
const amr::scalar_t& AmrMultiGridLevel<MatrixType, VectorType>::cellSize(lo_t dir) const {
frey_m's avatar
frey_m committed
195 196 197 198 199 200 201 202 203 204 205
    return dx_m[dir];
}


template <class MatrixType, class VectorType>
const amr::scalar_t* AmrMultiGridLevel<MatrixType, VectorType>::invCellSize() const {
    return invdx_m;
}


template <class MatrixType, class VectorType>
frey_m's avatar
frey_m committed
206
const amr::scalar_t& AmrMultiGridLevel<MatrixType, VectorType>::invCellSize(lo_t dir) const {
frey_m's avatar
frey_m committed
207 208 209 210
    return invdx_m[dir];
}


211 212 213 214 215 216 217
template <class MatrixType, class VectorType>
bool AmrMultiGridLevel<MatrixType, VectorType>::isValid(const AmrIntVect_t& iv) const {
    return ( iv.allGT(AmrIntVect_t(D_DECL(-1, -1, -1))) &&
             iv.allLT(AmrIntVect_t(D_DECL(nr_m[0], nr_m[1], nr_m[2]))) );
}


218
template <class MatrixType, class VectorType>
frey_m's avatar
frey_m committed
219 220
void AmrMultiGridLevel<MatrixType, VectorType>::buildMap(const Teuchos::RCP<comm_t>& comm,
                                                         const Teuchos::RCP<node_t>& node)
221 222
{
    
frey_m's avatar
frey_m committed
223
    go_t localNumElements = 0;
224
    
225
    Teuchos::Array<go_t> globalindices;
226 227 228 229 230 231 232 233 234 235 236 237 238
    
    for (amrex::MFIter mfi(grids, dmap, true); mfi.isValid(); ++mfi) {
        const amrex::Box&    tbx  = mfi.tilebox();
        const int* lo = tbx.loVect();
        const int* hi = tbx.hiVect();
        
        for (int i = lo[0]; i <= hi[0]; ++i) {
            for (int j = lo[1]; j <= hi[1]; ++j) {
#if AMREX_SPACEDIM == 3
                for (int k = lo[2]; k <= hi[2]; ++k) {
#endif
                    AmrIntVect_t iv(D_DECL(i, j, k));

frey_m's avatar
frey_m committed
239
                    go_t globalidx = serialize(iv);
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
                    
                    globalindices.push_back(globalidx);
                    
                    ++localNumElements;
#if AMREX_SPACEDIM == 3
                }
#endif
            }
        }
    }
    
    /*
     * create map that specifies which processor gets which data
     */
    
    // get smallest global index of this level
    amrex::Box bx = grids.minimalBox();
    const int* lo = bx.loVect();
    AmrIntVect_t lowcorner(D_DECL(lo[0], lo[1], lo[2]));
    
    // where to start indexing
261
    go_t baseIndex = serialize(lowcorner);
262 263
    
    // numGlobalElements == N
frey_m's avatar
frey_m committed
264
    go_t N = grids.numPts();
265 266 267
    
    map_p = Teuchos::rcp( new dmap_t(N, globalindices, baseIndex, comm, node) );
}