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
// 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/>.
//
21 22 23 24 25 26 27 28 29 30 31
#ifndef AMR_LAGRANGE_INTERPOLATER_H
#define AMR_LAGRANGE_INTERPOLATER_H

#include "AmrInterpolater.h"

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

32 33
#include "Ippl.h"

frey_m's avatar
frey_m committed
34 35
template <class Level>
class AmrLagrangeInterpolater : public AmrInterpolater<Level>
36 37 38
{
public:
    
39 40 41 42 43 44
    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;
45

46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
    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,
63 64 65
                 const basefab_t& fab,
                 umap_t& map,
                 const scalar_t& scale,
frey_m's avatar
frey_m committed
66
                 Level* mglevel);
67 68
    
    void coarse(const AmrIntVect_t& iv,
69 70 71
                umap_t& map,
                const scalar_t& scale,
                lo_t dir, lo_t shift, const basefab_t& rfab,
72
                const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
73
                Level* mglevel);
74 75
    
    void fine(const AmrIntVect_t& iv,
76 77 78
              umap_t& map,
              const scalar_t& scale,
              lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
79
              Level* mglevel);
80 81 82 83 84 85 86 87
    
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)
88 89 90
     * @param map with global matrix indices of fine level cells and
     * matrix entries of fine level cells (coefficients)
     * @param scale of matrix values
91 92 93 94 95 96 97 98
     * @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,
99 100 101
                      umap_t& map,
                      const scalar_t& scale,
                      lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
102
                      Level* mglevel);
103 104 105 106 107 108
    
    /*!
     * 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)
109 110 111
     * @param map with global matrix indices of fine level cells and
     * values matrix entries of fine level cells (coefficients)
     * @param scale of matrix values
112 113 114 115 116 117 118 119
     * @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,
120 121 122
                         umap_t& map,
                         const scalar_t& scale,
                         lo_t dir, lo_t shift,
frey_m's avatar
frey_m committed
123
                         Level* mglevel);
124 125 126 127
    
    /*!
     * First oder interpolation on coarse cell interface side
     * @param iv is the coarse cell at the interface (center cell of Laplacian)
128 129 130
     * @param map with global matrix indices of coarse level cells and
     * values matrix entries of coarse level cells (coefficients)
     * @param scale of matrix values
131 132 133 134 135 136 137 138 139 140
     * @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,
141 142 143
                      umap_t& map,
                      const scalar_t& scale,
                      lo_t dir, lo_t shift, const basefab_t& rfab,
144
                      const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
145
                      Level* mglevel);
146 147 148 149 150
    
    /*!
     * Second order interpolation on coarse cell interface side
     * 
     * @param iv is the coarse cell at the interface (center cell of Laplacian)
151 152 153
     * @param map with global matrix indices of coarse level cells and
     * values matrix entries of coarse level cells (coefficients)
     * @param scale of matrix values
154 155 156 157 158 159 160 161 162 163
     * @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,
164 165 166
                         umap_t& map,
                         const scalar_t& scale,
                         lo_t dir, lo_t shift, const basefab_t& rfab,
167
                         const AmrIntVect_t& riv,
frey_m's avatar
frey_m committed
168
                         Level* mglevel);
169 170
    
private:
171
#if AMREX_SPACEDIM == 3
172
    static constexpr qpattern_t qpattern_ms {
frey_m's avatar
frey_m committed
173 174 175 176 177 178 179 180 181
        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)
182 183 184 185 186 187 188 189 190
    };
    
    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
191 192

    // y_b   y_t
193 194 195 196
    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];
197 198 199 200 201 202 203 204 205 206
#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
207 208 209 210 211
};

#include "AmrLagrangeInterpolater.hpp"

#endif