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"
 
class AmrPartBunch;
 
// AMReX headers
#include <AMReX_AmrMesh.H>
#include <AMReX.H>
 
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