AmrLagrangeInterpolater.h 8.52 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4
//
// Class AmrLagrangeInterpolater
//   Lagrange interpolation for coarse-fine interfaces.
//
frey_m's avatar
frey_m committed
5
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
// All rights reserved
//
// 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/>.
//

22 23 24 25 26 27 28 29 30 31 32
#ifndef AMR_LAGRANGE_INTERPOLATER_H
#define AMR_LAGRANGE_INTERPOLATER_H

#include "AmrInterpolater.h"

#if AMREX_SPACEDIM == 3
    #include <bitset>
    #include <iterator>
    #include <array>
#endif

33 34
#include "Ippl.h"

frey_m's avatar
frey_m committed
35 36
template <class Level>
class AmrLagrangeInterpolater : public AmrInterpolater<Level>
37 38 39
{
public:
    
40 41 42 43 44 45
    typedef typename Level::go_t        go_t;
    typedef typename Level::lo_t        lo_t;
    typedef typename Level::scalar_t    scalar_t;
    typedef typename Level::umap_t      umap_t;
    typedef typename Level::basefab_t   basefab_t;
    typedef amr::AmrIntVect_t           AmrIntVect_t;
46

47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
    enum Order {
        LINEAR = 1,
        QUADRATIC
    };
    
#if AMREX_SPACEDIM == 3
    typedef std::bitset<25> qbits_t; ///< for checking the neighbour cells (quadratic)
    typedef std::bitset<9> lbits_t; ///< for checking the neighbour cells (linear)
    typedef std::array<unsigned int long, 9> qpattern_t;    ///< quadratic pattern
    typedef std::array<unsigned int long, 4> lpattern_t;    ///< linear pattern
#endif
    
public:
    
    AmrLagrangeInterpolater(Order order);
    
    void stencil(const AmrIntVect_t& iv,
64 65 66
                 const basefab_t& fab,
                 umap_t& map,
                 const scalar_t& scale,
frey_m's avatar
frey_m committed
67
                 Level* mglevel);
68 69
    
    void coarse(const AmrIntVect_t& iv,
70 71 72
                umap_t& map,
                const scalar_t& scale,
                lo_t dir, lo_t shift, const basefab_t& rfab,
73
                const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
74
                Level* mglevel);
75 76
    
    void fine(const AmrIntVect_t& iv,
77 78 79
              umap_t& map,
              const scalar_t& scale,
              lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
80
              Level* mglevel);
81 82 83 84 85 86 87 88
    
private:
    
    /*!
     * First order interpolation on fine cell interface side
     * 
     * @param iv is the fine ghost cell at the interface (on coarse cell that is not
     * refined)
89 90 91
     * @param map with global matrix indices of fine level cells and
     * matrix entries of fine level cells (coefficients)
     * @param scale of matrix values
92 93 94 95 96 97 98 99
     * @param dir direction of interface (0 "horizontal", 1 "vertical", 2 "longitudinal")
     * @param shift is either -1 or 1. If the refined coarse cell is on the left / lower / front
     * side, shift is equal to -1, otherwise the interface is on the right / upper / back side
     * and the value is 1.
     * @param mglevel used to get the global indices and refinement ratio among levels,
     * and boundary avlues at physical domain, e.g. Dirichlet, open BC
     */
    void fineLinear_m(const AmrIntVect_t& iv,
100 101 102
                      umap_t& map,
                      const scalar_t& scale,
                      lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
103
                      Level* mglevel);
104 105 106 107 108 109
    
    /*!
     * Second order interpolation on fine cell interface side
     * 
     * @param iv is the fine ghost cell at the interface (on coarse cell that is not
     * refined)
110 111 112
     * @param map with global matrix indices of fine level cells and
     * values matrix entries of fine level cells (coefficients)
     * @param scale of matrix values
113 114 115 116 117 118 119 120
     * @param dir direction of interface (0 "horizontal", 1 "vertical", 2 "longitudinal")
     * @param shift is either -1 or 1. If the refined coarse cell is on the left / lower / front
     * side, shift is equal to -1, otherwise the interface is on the right / upper / back side
     * and the value is 1.
     * @param mglevel used to get the global indices and refinement ratio among levels,
     * and boundary avlues at physical domain, e.g. Dirichlet, open BC
     */
    void fineQuadratic_m(const AmrIntVect_t& iv,
121 122 123
                         umap_t& map,
                         const scalar_t& scale,
                         lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
124
                         Level* mglevel);
125 126 127 128
    
    /*!
     * First oder interpolation on coarse cell interface side
     * @param iv is the coarse cell at the interface (center cell of Laplacian)
129 130 131
     * @param map with global matrix indices of coarse level cells and
     * values matrix entries of coarse level cells (coefficients)
     * @param scale of matrix values
132 133 134 135 136 137 138 139 140 141
     * @param dir direction of interface (0 "horizontal", 1 "vertical", 2 "longitudinal")
     * @param shift is either -1 or 1. If the refined coarse cell is on the left / lower / front
     * side, shift is equal to -1, otherwise the interface is on the right / upper / back side
     * and the value is 1.
     * @param ba contains all coarse cells that got refined
     * @param riv is the fine cell at the interface
     * @param mglevel used to get the global indices and refinement ratio among levels,
     * and boundary values at physical domain, e.g. Dirichlet, open BC
     */
    void crseLinear_m(const AmrIntVect_t& iv,
142 143 144
                      umap_t& map,
                      const scalar_t& scale,
                      lo_t dir, lo_t shift, const basefab_t& rfab,
145
                      const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
146
                      Level* mglevel);
147 148 149 150 151
    
    /*!
     * Second order interpolation on coarse cell interface side
     * 
     * @param iv is the coarse cell at the interface (center cell of Laplacian)
152 153 154
     * @param map with global matrix indices of coarse level cells and
     * values matrix entries of coarse level cells (coefficients)
     * @param scale of matrix values
155 156 157 158 159 160 161 162 163 164
     * @param dir direction of interface (0 "horizontal", 1 "vertical", 2 "longitudinal")
     * @param shift is either -1 or 1. If the refined coarse cell is on the left / lower / front
     * side, shift is equal to -1, otherwise the interface is on the right / upper / back side
     * and the value is 1.
     * @param ba contains all coarse cells that got refined
     * @param riv is the fine cell at the interface
     * @param mglevel used to get the global indices and refinement ratio among levels,
     * and boundary values at physical domain, e.g. Dirichlet, open BC
     */
    void crseQuadratic_m(const AmrIntVect_t& iv,
165 166 167
                         umap_t& map,
                         const scalar_t& scale,
                         lo_t dir, lo_t shift, const basefab_t& rfab,
168
                         const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
169
                         Level* mglevel);
170 171
    
private:
172
#if AMREX_SPACEDIM == 3
173
    static constexpr qpattern_t qpattern_ms {
frey_m's avatar
frey_m committed
174 175 176 177 178 179 180 181 182
        473536,                             ///< cross pattern (case 0)
        14798,                              ///< T pattern (case 1)
        236768,                             ///< right hammer pattern (case 2)
        15153152,                           ///< T on head pattern (case 3)
        947072,                             ///< left hammer pattern (case 4)
        29596,                              ///< upper left corner pattern (case 5)
        7399,                               ///< upper right corner pattern (case 6)
        7576576,                            ///< mirrored L pattern (case 7)
        30306304                            ///< L pattern (case 8)
183 184 185 186 187 188 189 190 191
    };
    
    static constexpr lpattern_t lpattern_ms {
        27,                                 ///< corner top right pattern
        216,                                ///< corner bottom right pattern
        432,                                ///< corner bottom left pattern
        54                                  ///< corner top left pattern
    };
#endif
192 193

    // y_b   y_t
194 195 196 197
    static const scalar_t lookup1a_ms[2];
    static const scalar_t lookup2a_ms[2];
    static const scalar_t lookup1b_ms[2];
    static const scalar_t lookup2b_ms[2];
198 199 200 201 202 203 204 205 206 207
#if AMREX_SPACEDIM == 3
    static const scalar_t lookup3_ms[2];
    static const scalar_t lookup3r_ms[2];
    static const scalar_t lookup4_ms[2];
    static const scalar_t lookup4r_ms[2];
    static const scalar_t lookup5_ms[2];
    static const scalar_t lookup5r_ms[2];
    static const scalar_t lookup6_ms;
    static const scalar_t factor_ms;
#endif
208 209 210 211 212
};

#include "AmrLagrangeInterpolater.hpp"

#endif