AmrBoxLib.h

Go to the documentation of this file.

//
// Class AmrBoxLib
//   Concrete AMR object. It is based on the AMReX library
//   (cf. https://amrex-codes.github.io/ or https://ccse.lbl.gov/AMReX/).
//   AMReX is the successor of BoxLib. This class represents the interface
//   to AMReX and the AMR framework in OPAL. It implements the functions of
//   the AmrObject class.
//
// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// 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/>.
//
#ifndef AMR_BOXLIB_H
#define AMR_BOXLIB_H

#include "Amr/AmrObject.h"

// AMReX headers
#include <AMReX_AmrMesh.H>
#include <AMReX.H>

#include <map>
#include <vector>

class AmrPartBunch;

class AmrBoxLib: public AmrObject,
                 public amrex::AmrMesh {

public:
    typedef amr::AmrField_t                 AmrField_t;
    typedef amr::AmrScalarFieldContainer_t  AmrScalarFieldContainer_t;
    typedef amr::AmrVectorFieldContainer_t  AmrVectorFieldContainer_t;
    typedef amr::AmrGeomContainer_t         AmrGeomContainer_t;
    typedef amr::AmrGridContainer_t         AmrGridContainer_t;
    typedef amr::AmrProcMapContainer_t      AmrProcMapContainer_t;
    typedef amr::AmrDomain_t                AmrDomain_t;
    typedef amr::AmrIntArray_t              AmrIntArray_t;
    typedef amr::AmrReal_t                  AmrReal_t;
    typedef amr::AmrGrid_t                  AmrGrid_t;
    typedef amr::AmrProcMap_t               AmrProcMap_t;
    typedef amr::AmrGeometry_t              AmrGeometry_t;
    typedef amr::AmrIntVect_t               AmrIntVect_t;

    typedef amrex::FArrayBox                FArrayBox_t;
    typedef amrex::Box                      Box_t;
    typedef amrex::TagBox                   TagBox_t;
    typedef amrex::TagBoxArray              TagBoxArray_t;
    typedef amrex::MFIter                   MFIter_t;

    enum Strategy {
        RANK_ZERO = 0,      // all grids to processor zero
        PFC       = 1,
        RANDOM    = 2,
        KNAPSACK  = 3
    };

public:
    AmrBoxLib(const AmrDomain_t& domain,
              const AmrIntArray_t& nGridPts,
              int maxLevel,
              AmrPartBunch* bunch_p);

    static std::unique_ptr<AmrBoxLib> create(const AmrInfo& info,
                                             AmrPartBunch* bunch_p);

    void regrid(double time);

    void getGridStatistics(std::map<int, long>& gridPtsPerCore,
                           std::vector<int>& gridsPerLevel) const;

    void initFineLevels();

    VectorPair_t getEExtrema();

    double getRho(int x, int y, int z);

    void computeSelfFields();

    void computeSelfFields(int bin);

    void computeSelfFields_cycl(double gamma);

    void computeSelfFields_cycl(int bin);

    void updateMesh();

    const Vector_t& getMeshScaling() const;

    Vektor<int, 3> getBaseLevelGridPoints() const;

    const int& maxLevel() const;
    const int& finestLevel() const;

    double getT() const;

    void redistributeGrids(int how);


protected:
    /*
     * AmrMesh functions
     */

    void RemakeLevel (int lev, AmrReal_t time,
                      const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap);

    void MakeNewLevel (int lev, AmrReal_t time,
                       const AmrGrid_t& new_grids, const AmrProcMap_t& new_dmap);

    void ClearLevel(int lev);

    virtual void ErrorEst(int lev, TagBoxArray_t& tags,
                          AmrReal_t time, int ngrow) override;

    void MakeNewLevelFromScratch(int lev, AmrReal_t time,
                                 const AmrGrid_t& ba,
                                 const AmrProcMap_t& dm);

    void MakeNewLevelFromCoarse(int lev, AmrReal_t time,
                                const AmrGrid_t& ba,
                                const AmrProcMap_t& dm);


private:
    void doRegrid_m(int lbase, double time);

    void preRegrid_m();

    void postRegrid_m(int old_finest);

    double solvePoisson_m();

    /* ATTENTION
     * The tagging routines assume the particles to be in the
     * AMR domain, i.e. [-1, 1]^3
     */

    void tagForChargeDensity_m(int lev, TagBoxArray_t& tags,
                               AmrReal_t time, int ngrow);

    void tagForPotentialStrength_m(int lev, TagBoxArray_t& tags,
                                   AmrReal_t time, int ngrow);

    void tagForEfield_m(int lev, TagBoxArray_t& tags,
                        AmrReal_t time, int ngrow);

    void tagForMomenta_m(int lev, TagBoxArray_t& tags,
                         AmrReal_t time, int ngrow);

    void tagForMaxNumParticles_m(int lev, TagBoxArray_t& tags,
                                 AmrReal_t time, int ngrow);

    void tagForMinNumParticles_m(int lev, TagBoxArray_t& tags,
                                 AmrReal_t time, int ngrow);

    void initBaseLevel_m(const AmrIntArray_t& nGridPts);

    static void initParmParse_m(const AmrInfo& info, AmrLayout_t* layout_p);

    void fillPhysbc_m(AmrField_t& mf, int lev = 0);

//     void gradient(int lev) {
//         
//         phi_m[lev]->FillBoundary(geom[lev].periodicity());
//         
//         for (MFIter_t mfi(*phi_m[lev], true); mfi.isValid(); ++mfi) {
//             const Box_t&  tilebx  = mfi.tilebox();
//             FArrayBox_t&  fab = (*phi_m[lev])[mfi];
//             FArrayBox_t& efab = (*efield_m[lev])[mfi];
//             
//             const int*  tlo     = tilebx.loVect();
//             const int*  thi     = tilebx.hiVect();
//             
//             for (int i = tlo[0]; i <= thi[0]; ++i) {
//                 for (int j = tlo[1]; j <= thi[1]; ++j) {
//                     for (int k = tlo[2]; k <= thi[2]; ++k) {
//                         
//                         amrex::IntVect iv(D_DECL(i,j,k));
//                         
//                         // x-field
//                         amrex::IntVect liv(D_DECL(i-1,j,k));
//                         amrex::IntVect riv(D_DECL(i+1,j,k));
//                         efab(iv, 0) = 0.5 * (fab(liv) - fab(riv));
//                         
//                         // y-field
//                         liv = amrex::IntVect(D_DECL(i,j-1,k));
//                         riv = amrex::IntVect(D_DECL(i,j+1,k));
//                         efab(iv, 1) = 0.5 * (fab(liv) - fab(riv));
//                         
//                         // z-field
//                         liv = amrex::IntVect(D_DECL(i,j-1,k));
//                         riv = amrex::IntVect(D_DECL(i,j+1,k));
//                         efab(iv, 2) = 0.5 * (fab(liv) - fab(riv));
//                     }
//                 }
//             }
//         }
//     }


private:
    AmrPartBunch* bunch_mp;

    // the layout of the bunch
    AmrLayout_t* layout_mp;

    AmrScalarFieldContainer_t rho_m;

    AmrScalarFieldContainer_t phi_m;

    AmrVectorFieldContainer_t efield_m;

    Vector_t meshScaling_m;

    // only used in case of potential and efield tagging
    std::vector<bool> isFirstTagging_m;

    bool isPoissonSolved_m;
};

#endif