AmrOpenBoundary.h 11 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4
//
// Class AmrOpenBoundary
//   Open boundary condition
//
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
#ifndef AMR_OPEN_BOUNDARY_H
#define AMR_OPEN_BOUNDARY_H

#include "AmrBoundary.h"

frey_m's avatar
frey_m committed
27 28
#include "Utilities/OpalException.h"

frey_m's avatar
frey_m committed
29 30
template <class Level>
class AmrOpenBoundary : public AmrBoundary<Level> {
31 32

public:
frey_m's avatar
frey_m committed
33 34
    typedef typename Level::umap_t      umap_t;
    typedef typename Level::lo_t        lo_t;
frey_m's avatar
frey_m committed
35
    typedef typename Level::go_t        go_t;
frey_m's avatar
frey_m committed
36 37
    typedef typename Level::scalar_t    scalar_t;
    typedef amr::AmrIntVect_t           AmrIntVect_t;
38 39 40
    
public:
    
frey_m's avatar
frey_m committed
41 42 43 44 45
    AmrOpenBoundary()
        : AmrBoundary<Level>(4)
        , order_m(ABC::Robin)
        , dist_m(1.7)
    { }
46 47
    
    void apply(const AmrIntVect_t& iv,
48 49 50
               const lo_t& dir,
               umap_t& map,
               const scalar_t& value,
frey_m's avatar
frey_m committed
51
               Level* mglevel,
frey_m's avatar
frey_m committed
52
               const go_t* nr);
53
    
frey_m's avatar
frey_m committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 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 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
    enum ABC {
        Zeroth,
        First,
        Second,
        Third,
        Robin
    };
    
private:
    int order_m;
    double dist_m;
    
    scalar_t coordinate_m(const AmrIntVect_t& iv,
                          const lo_t& dir,
                          Level* mglevel,
                          const go_t* nr);
    
    /*!
     * Robin boundary condition
     */
    void robin_m(const AmrIntVect_t& iv,
                 const lo_t& dir,
                 umap_t& map,
                 const scalar_t& value,
                 Level* mglevel,
                 const go_t* nr);
    
    /*!
     * Asymptotic boundary condition 0th order (ABC0)
     */
    void abc0_m(const AmrIntVect_t& iv,
                const lo_t& dir,
                umap_t& map,
                const scalar_t& value,
                Level* mglevel,
                const go_t* nr);
    
    
    /*!
     * Asymptotic boundary condition 1st order (ABC1)
     */
    void abc1_m(const AmrIntVect_t& iv,
                const lo_t& dir,
                umap_t& map,
                const scalar_t& value,
                Level* mglevel,
                const go_t* nr);
    
    /*!
     * Asymptotic boundary condition 2nd order (ABC2)
     */
    void abc2_m(const AmrIntVect_t& iv,
                const lo_t& dir,
                umap_t& map,
                const scalar_t& value,
                Level* mglevel,
                const go_t* nr);
    
    
    /*!
     * Asymptotic boundary condition 3rd order (ABC3)
     */
    void abc3_m(const AmrIntVect_t& iv,
                const lo_t& dir,
                umap_t& map,
                const scalar_t& value,
                Level* mglevel,
                const go_t* nr);
    
123 124 125
};


frey_m's avatar
frey_m committed
126 127 128 129 130 131
template <class Level>
void AmrOpenBoundary<Level>::apply(const AmrIntVect_t& iv,
                                   const lo_t& dir,
                                   umap_t& map,
                                   const scalar_t& value,
                                   Level* mglevel,
frey_m's avatar
frey_m committed
132
                                   const go_t* nr)
133
{
frey_m's avatar
frey_m committed
134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166
    switch ( order_m ) {
        case ABC::Third:
        {
            this->abc3_m(iv, dir, map, value, mglevel, nr);
            break;
        }
        case ABC::Second:
        {
            this->abc2_m(iv, dir, map, value, mglevel, nr);
            break;
        }
        case ABC::First:
        {
            this->abc1_m(iv, dir, map, value, mglevel, nr);
            break;
        }
        case ABC::Zeroth:
        {
            this->abc0_m(iv, dir, map, value, mglevel, nr);
            break;
        }
        case ABC::Robin:
        {
            this->robin_m(iv, dir, map, value, mglevel, nr);
            break;
        }
        default:
        {
            throw OpalException("AmrOpenBoundary<Level>::apply()",
                                "Order not supported.");
        }
    }
}
167

frey_m's avatar
frey_m committed
168 169 170 171 172 173 174 175 176

template <class Level>
void AmrOpenBoundary<Level>::robin_m(const AmrIntVect_t& iv,
                                     const lo_t& dir,
                                     umap_t& map,
                                     const scalar_t& value,
                                     Level* mglevel,
                                     const go_t* nr)
{
177
    AmrIntVect_t niv = iv;
178
    AmrIntVect_t n2iv = iv;
179
    
frey_m's avatar
frey_m committed
180 181 182
    if ( iv[dir] == -1 ) {
        // lower boundary
        niv[dir] = 0;
183 184
        n2iv[dir] = 1;
    } else {
frey_m's avatar
frey_m committed
185 186
        // upper boundary        
        niv[dir] = nr[dir] - 1;
187
        n2iv[dir] = nr[dir] - 2;
188 189
    }
    
frey_m's avatar
frey_m committed
190
    // correct cell size
frey_m's avatar
frey_m committed
191
    scalar_t h = mglevel->cellSize(dir);
frey_m's avatar
frey_m committed
192 193 194 195 196 197 198
    
    // artificial distance
    scalar_t d = dist_m;;
    
    map[mglevel->serialize(niv)] -= 2.0 * h / d * value;
    map[mglevel->serialize(n2iv)] += value;
}
frey_m's avatar
frey_m committed
199

frey_m's avatar
frey_m committed
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267

template <class Level>
void AmrOpenBoundary<Level>::abc0_m(const AmrIntVect_t& iv,
                                    const lo_t& dir,
                                    umap_t& map,
                                    const scalar_t& value,
                                    Level* mglevel,
                                    const go_t* nr)
{
    // ABC0 == Dirichlet BC
    AmrIntVect_t niv = iv;
    
    if ( iv[dir] == -1 ) {
        // lower boundary
        niv[dir] = 0;
        
    } else {
        // upper boundary        
        niv[dir] = nr[dir] - 1;
    }
    
    map[mglevel->serialize(niv)] -= value;
}


template <class Level>
void AmrOpenBoundary<Level>::abc1_m(const AmrIntVect_t& iv,
                                    const lo_t& dir,
                                    umap_t& map,
                                    const scalar_t& value,
                                    Level* mglevel,
                                    const go_t* nr)
{
    // ABC1 == Robin BC (with normal derivatives) (i.e. Dirichlet BC + Neumann BC)
    AmrIntVect_t niv = iv;
    AmrIntVect_t n2iv = iv;
    
    scalar_t sign = 1.0;
    
    if ( iv[dir] == -1 ) {
        // lower boundary
        niv[dir] = 0;
        n2iv[dir] = 1;
        
        sign = -1.0;
        
    } else {
        // upper boundary        
        niv[dir] = nr[dir] - 1;
        n2iv[dir] = nr[dir] - 2;
    }
    
    // correct cell size
    scalar_t h = mglevel->cellSize(dir);
    
    // coordinate of inner cell
    scalar_t x = this->coordinate_m(niv, 0, mglevel, nr) + sign * dist_m * h;
    scalar_t y = this->coordinate_m(niv, 1, mglevel, nr) + sign * dist_m * h;
#if AMREX_SPACEDIM == 3
    scalar_t z = this->coordinate_m(niv, 2, mglevel, nr) + sign * dist_m * h;
#endif
    scalar_t r = std::sqrt( AMREX_D_TERM(x * x, + y * y, + z * z) );
    
//     // artificial distance
//     scalar_t d = dist_m * nr[dir] * h; // + r;
    scalar_t d = r;
    
    map[mglevel->serialize(niv)] -= 2.0 * h / d * value;
frey_m's avatar
frey_m committed
268
    map[mglevel->serialize(n2iv)] += value;
269 270
}

frey_m's avatar
frey_m committed
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397

template <class Level>
void AmrOpenBoundary<Level>::abc2_m(const AmrIntVect_t& iv,
                                    const lo_t& dir,
                                    umap_t& map,
                                    const scalar_t& value,
                                    Level* mglevel,
                                    const go_t* nr)
{
    AmrIntVect_t niv = iv;
    AmrIntVect_t n2iv = iv;
    
    scalar_t sign = 1.0;
    
    if ( iv[dir] == -1 ) {
        // lower boundary
        niv[dir] = 0;
        n2iv[dir] = 1;
        
        sign = -1.0;
        
    } else {
        // upper boundary        
        niv[dir] = nr[dir] - 1;
        n2iv[dir] = nr[dir] - 2;
    }
    
    // correct cell size
    scalar_t h = mglevel->cellSize(dir);
    
    // coordinate of inner cell
    scalar_t x = this->coordinate_m(niv, 0, mglevel, nr) + sign * dist_m * h;
    scalar_t y = this->coordinate_m(niv, 1, mglevel, nr) + sign * dist_m * h;
#if AMREX_SPACEDIM == 3
    scalar_t z = this->coordinate_m(niv, 2, mglevel, nr) + sign * dist_m * h;
#endif
    scalar_t r = std::sqrt( AMREX_D_TERM(x * x, + y * y, + z * z) );
    
//     scalar_t d = dist_m * nr[dir] * h; // + r;
    scalar_t d = r;
    
    map[mglevel->serialize(niv)] += 2.0 * (d * d - h * h) / ( d * d + 2.0 * h * d ) * value;
    map[mglevel->serialize(n2iv)] += (2.0 * h - d) / (2.0 * h + d) * value;
}


template <class Level>
void AmrOpenBoundary<Level>::abc3_m(const AmrIntVect_t& iv,
                                    const lo_t& dir,
                                    umap_t& map,
                                    const scalar_t& value,
                                    Level* mglevel,
                                    const go_t* nr)
{
    AmrIntVect_t niv = iv;
    AmrIntVect_t n1iv = iv;
    AmrIntVect_t n2iv = iv;
    AmrIntVect_t n3iv = iv;
    AmrIntVect_t n4iv = iv;
    
    scalar_t sign = 1.0;
    
    if ( iv[dir] == -1 ) {
        // lower boundary
        niv[dir] = 0;
        n1iv[dir] = 1;
        n2iv[dir] = 2;
        n3iv[dir] = 3;
        n4iv[dir] = 4;
        
        
        
    } else {
        // upper boundary        
        niv[dir] = nr[dir] - 1;
        n1iv[dir] = nr[dir] - 2;
        n2iv[dir] = nr[dir] - 3;
        n3iv[dir] = nr[dir] - 4;
        n4iv[dir] = nr[dir] - 5;
        
        sign *= -1.0;
    }
    
    scalar_t h = mglevel->cellSize(dir);
    
//     scalar_t d = dist_m * nr[dir] * h; // + r;
    
    // coordinate of inner cell
    scalar_t x = this->coordinate_m(niv, 0, mglevel, nr);
    scalar_t y = this->coordinate_m(niv, 1, mglevel, nr);
#if AMREX_SPACEDIM == 3
    scalar_t z = this->coordinate_m(niv, 2, mglevel, nr);
#endif
    
    scalar_t d = std::sqrt( AMREX_D_TERM(x * x, + y * y, + z * z) );
    
    scalar_t hd = h + d;
    
    map[mglevel->serialize(niv)]  += (2.0 * d  / hd - 2.0 * h * h / (3.0 * d * hd) + sign * 5.0 * d * d / (18.0 * h * hd)) * value;
    map[mglevel->serialize(n1iv)] += (h / hd - d / hd - sign * d * d / ( hd * h) ) * value;
    map[mglevel->serialize(n2iv)] += sign * 4.0 * d * d / (3.0 * h * hd) * value;
    map[mglevel->serialize(n3iv)] -= sign * 7.0 * d * d / (9.0 * h * hd) * value;
    map[mglevel->serialize(n4iv)] += sign * d * d / (6.0 * h * hd) * value;
}


template <class Level>
typename AmrOpenBoundary<Level>::scalar_t
AmrOpenBoundary<Level>::coordinate_m(const AmrIntVect_t& iv,
                                     const lo_t& dir,
                                     Level* mglevel,
                                     const go_t* nr)
{
    /*
     * subtract / add half cell length in order to get
     * cell centered position
     */
    scalar_t h = mglevel->cellSize(dir);
    
    scalar_t lower = mglevel->geom.ProbLo(dir) + 0.5 * h;
    scalar_t upper = mglevel->geom.ProbHi(dir) - 0.5 * h;
    
    scalar_t m = (upper - lower) / ( nr[dir] - 1 );
    
    return m * iv[dir] + lower;
}

398
#endif