FMGPoissonSolver.cpp 11.8 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
//
// Class FMGPoissonSolver
//   Interface to the Fortran based AMR Poisson multigrid solver of BoxLib.
//   The Fortran solver is part of AMReX till version 18.07
//   (Link: https://github.com/AMReX-Codes/amrex/archive/18.07.tar.gz)
//
// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, 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/>.
//
frey_m's avatar
frey_m committed
23 24
#include "FMGPoissonSolver.h"

25 26
#include "Utility/PAssert.h"

27 28
#include "Utilities/OpalException.h"

frey_m's avatar
frey_m committed
29
#include <AMReX_ParmParse.H>
30
#include <AMReX_Interpolater.H>
frey_m's avatar
frey_m committed
31

frey_m's avatar
frey_m committed
32 33
FMGPoissonSolver::FMGPoissonSolver(AmrBoxLib* itsAmrObject_p)
    : AmrPoissonSolver<AmrBoxLib>(itsAmrObject_p),
frey_m's avatar
frey_m committed
34 35
      reltol_m(1.0e-15),
      abstol_m(1.0e-9)
frey_m's avatar
frey_m committed
36 37
{
    // Dirichlet boundary conditions are default
38
    for (int d = 0; d < AMREX_SPACEDIM; ++d) {
frey_m's avatar
frey_m committed
39 40 41
        bc_m[2 * d]     = MGT_BC_DIR;
        bc_m[2 * d + 1] = MGT_BC_DIR;
    }
kraus's avatar
kraus committed
42

frey_m's avatar
frey_m committed
43
    this->initParameters_m();
frey_m's avatar
frey_m committed
44 45
}

46 47 48
void FMGPoissonSolver::solve(AmrScalarFieldContainer_t& rho,
                             AmrScalarFieldContainer_t& phi,
                             AmrVectorFieldContainer_t& efield,
frey_m's avatar
frey_m committed
49
                             unsigned short baseLevel,
frey_m's avatar
frey_m committed
50
                             unsigned short finestLevel,
frey_m's avatar
frey_m committed
51
                             bool /*prevAsGuess*/)
frey_m's avatar
frey_m committed
52
{
frey_m's avatar
frey_m committed
53
    const GeomContainer_t& geom = itsAmrObject_mp->Geom();
54

frey_m's avatar
frey_m committed
55
    if (AmrGeometry_t::isAllPeriodic()) {
56
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
frey_m's avatar
frey_m committed
57 58 59
            bc_m[2 * d]     = MGT_BC_PER;
            bc_m[2 * d + 1] = MGT_BC_PER;
        }
frey_m's avatar
frey_m committed
60
    } else if ( AmrGeometry_t::isAnyPeriodic() ) {
61
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
frey_m's avatar
frey_m committed
62
            if ( AmrGeometry_t::isPeriodic(d) ) {
frey_m's avatar
frey_m committed
63 64 65 66 67
                bc_m[2 * d]     = MGT_BC_PER;
                bc_m[2 * d + 1] = MGT_BC_PER;
            }
        }
    }
kraus's avatar
kraus committed
68

69
    amrex::Vector< AmrScalarFieldContainer_t > grad_phi_edge(rho.size());
kraus's avatar
kraus committed
70

71 72
    amrex::Vector< AmrField_t > tmp(rho.size());

73 74
    PAssert(baseLevel <= finestLevel);
    PAssert(finestLevel < (int)rho.size());
75

frey_m's avatar
frey_m committed
76
    for (int lev = baseLevel; lev <= finestLevel ; ++lev) {
frey_m's avatar
frey_m committed
77
        const AmrProcMap_t& dmap = rho[lev]->DistributionMap();
kraus's avatar
kraus committed
78
        grad_phi_edge[lev].resize(AMREX_SPACEDIM);
79
        tmp[lev].define(rho[lev]->boxArray(), dmap, AMREX_SPACEDIM, 1);
80

81
        for (int n = 0; n < AMREX_SPACEDIM; ++n) {
kraus's avatar
kraus committed
82
            AmrGrid_t ba = rho[lev]->boxArray();
frey_m's avatar
frey_m committed
83
            grad_phi_edge[lev][n].reset(new AmrField_t(ba.surroundingNodes(n), dmap, 1, 1));
kraus's avatar
kraus committed
84
            grad_phi_edge[lev][n]->setVal(0.0, 1);
frey_m's avatar
frey_m committed
85 86
        }
    }
87

frey_m's avatar
frey_m committed
88
    for (int i = 0; i <= finestLevel; ++i) {
frey_m's avatar
frey_m committed
89
        phi[i]->setVal(0.0, 1);
90
        tmp[i].setVal(0.0, 1);
frey_m's avatar
frey_m committed
91
    }
kraus's avatar
kraus committed
92

frey_m's avatar
frey_m committed
93 94 95
    double residNorm = this->solveWithF90_m(amrex::GetVecOfPtrs(rho),
                                            amrex::GetVecOfPtrs(phi),
                                            amrex::GetVecOfVecOfPtrs(grad_phi_edge),
96 97 98
                                            geom,
                                            baseLevel,
                                            finestLevel);
99

frey_m's avatar
frey_m committed
100
    if ( residNorm > abstol_m ) {
101 102
        std::stringstream ss;
        ss << "Residual norm: " << std::setprecision(16) << residNorm
frey_m's avatar
frey_m committed
103
           << " > " << abstol_m << " (absolute tolerance)";
104
        throw OpalException("FMGPoissonSolver::solve()",
105 106
                            "Multigrid solver did not converge. " + ss.str());
    }
kraus's avatar
kraus committed
107

frey_m's avatar
frey_m committed
108
    for (int lev = baseLevel; lev <= finestLevel; ++lev) {
109 110 111
        amrex::average_face_to_cellcenter(tmp[lev],
                                          amrex::GetVecOfConstPtrs(grad_phi_edge[lev]),
                                          geom[lev]);
kraus's avatar
kraus committed
112

113 114 115 116
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            efield[lev][j]->setVal(0.0, 1);
            amr::AmrField_t::Copy(*(efield[lev][j].get()),
                                  tmp[lev], j, 0, 1, 1);
117 118 119 120
            efield[lev][j]->FillBoundary(0, 1, geom[lev].periodicity());
            // we need also minus sign due to \vec{E} = - \nabla\phi
            efield[lev][j]->mult(-1.0, 0, 1);
        }
frey_m's avatar
frey_m committed
121 122 123 124
    }
}


125
double FMGPoissonSolver::getXRangeMin(unsigned short level) {
frey_m's avatar
frey_m committed
126
    return itsAmrObject_mp->Geom(level).ProbLo(0);
frey_m's avatar
frey_m committed
127 128 129
}


130
double FMGPoissonSolver::getXRangeMax(unsigned short level) {
frey_m's avatar
frey_m committed
131
    return itsAmrObject_mp->Geom(level).ProbHi(0);
frey_m's avatar
frey_m committed
132 133 134
}


135
double FMGPoissonSolver::getYRangeMin(unsigned short level) {
frey_m's avatar
frey_m committed
136
    return itsAmrObject_mp->Geom(level).ProbLo(1);
frey_m's avatar
frey_m committed
137 138 139
}


140
double FMGPoissonSolver::getYRangeMax(unsigned short level) {
frey_m's avatar
frey_m committed
141
    return itsAmrObject_mp->Geom(level).ProbHi(1);
frey_m's avatar
frey_m committed
142 143 144
}


145
double FMGPoissonSolver::getZRangeMin(unsigned short level) {
frey_m's avatar
frey_m committed
146
    return itsAmrObject_mp->Geom(level).ProbLo(2);
frey_m's avatar
frey_m committed
147 148 149
}


150
double FMGPoissonSolver::getZRangeMax(unsigned short level) {
frey_m's avatar
frey_m committed
151
    return itsAmrObject_mp->Geom(level).ProbHi(2);
frey_m's avatar
frey_m committed
152 153 154 155 156
}


Inform &FMGPoissonSolver::print(Inform &os) const {
    os << "* ************* F M G P o i s s o n S o l v e r ************************************ " << endl
157
       << "* relative tolerance " << reltol_m << '\n'
frey_m's avatar
frey_m committed
158 159 160 161 162 163
       << "* absolute tolerance " << abstol_m << '\n' << endl
       << "* ******************************************************************** " << endl;
    return os;
}


frey_m's avatar
frey_m committed
164 165 166 167
void FMGPoissonSolver::initParameters_m() {
    /* The information about the parameters is copied from the
     * BoxLib documentation chapter 5 on linear solvers and from
     * the code itself.
kraus's avatar
kraus committed
168
     *
frey_m's avatar
frey_m committed
169 170 171 172 173 174
     * Some paramters that have to be given to the solver using
     * AMReX_ParmParse.
     * "doc" tells how the variable is called in
     *  - amrex/Src/LinearSolvers/F_MG/cc_mg_cpp.f90
     *  - amrex/Src/LinearSolvers/F_MG/mg_tower.f90 (conains defaults)
     *  - amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
kraus's avatar
kraus committed
175
     *
frey_m's avatar
frey_m committed
176
     */
kraus's avatar
kraus committed
177

frey_m's avatar
frey_m committed
178
    amrex::ParmParse pp_mg("mg");
kraus's avatar
kraus committed
179

frey_m's avatar
frey_m committed
180 181
    // maximum number of multigrid cycles (doc: max_iter)
    pp_mg.add("maxiter", 200);
kraus's avatar
kraus committed
182

frey_m's avatar
frey_m committed
183 184
    // maximum number of iterations for the bottom solver (doc: bottom_max_iter)
    pp_mg.add("maxiter_b", 200);
kraus's avatar
kraus committed
185

frey_m's avatar
frey_m committed
186 187
    // number of smoothings at each level on the way DOWN the V-cycle (doc: nu1)
    pp_mg.add("nu_1", 2);
kraus's avatar
kraus committed
188

frey_m's avatar
frey_m committed
189 190
    // number of smoothings at each level on the way UP the V-cycle (doc: nu2)
    pp_mg.add("nu_2", 2);
kraus's avatar
kraus committed
191

frey_m's avatar
frey_m committed
192 193
    // number of smoothing before and after the bottom solver (doc: nub)
    pp_mg.add("nu_b", 0);
kraus's avatar
kraus committed
194

frey_m's avatar
frey_m committed
195 196
    // number of smoothing ... ?
    pp_mg.add("nu_f", 8);
kraus's avatar
kraus committed
197

frey_m's avatar
frey_m committed
198
    // verbosity of the multigrid solver. Higher numbers give more verbosity (doc: verbose)
199
    pp_mg.add("v"   , 0);
kraus's avatar
kraus committed
200 201


frey_m's avatar
frey_m committed
202 203
    // see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
    pp_mg.add("usecg", 1);
kraus's avatar
kraus committed
204

frey_m's avatar
frey_m committed
205 206
    // bottom_solver_eps, see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
    pp_mg.add("rtol_b", 0.0001);
kraus's avatar
kraus committed
207

frey_m's avatar
frey_m committed
208 209
    // see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
    pp_mg.add("cg_solver", 1);
kraus's avatar
kraus committed
210 211

    // maximum number of levels (max_nlevel)
frey_m's avatar
frey_m committed
212
    pp_mg.add("numLevelsMAX", 1024);
kraus's avatar
kraus committed
213

frey_m's avatar
frey_m committed
214 215 216 217 218 219 220
    /* MG_SMOOTHER_GS_RB  = 1   (red-black Gauss-Seidel)
     * MG_SMOOTHER_JACOBI = 2   (Jacobi)
     * MG_SMOOTHER_MINION_CROSS = 5
     * MG_SMOOTHER_MINION_FULL = 6
     * MG_SMOOTHER_EFF_RB = 7
     */
    pp_mg.add("smoother", 1);
kraus's avatar
kraus committed
221

frey_m's avatar
frey_m committed
222
    /* The type of multigrid to do
kraus's avatar
kraus committed
223
     *
frey_m's avatar
frey_m committed
224 225 226 227 228 229
     * MG_FCycle = 1    (F, full multigrid)
     * MG_WCycle = 2    (W-cycles)
     * MG_VCycle = 3    (V-cycles)
     * MG_FVCycle = 4   (F + V cycles)
     */
    pp_mg.add("cycle_type", 1);
kraus's avatar
kraus committed
230

frey_m's avatar
frey_m committed
231
    /* the type of bottom solver to use
kraus's avatar
kraus committed
232
     *
frey_m's avatar
frey_m committed
233 234 235 236
     * BiCG:    biconjugate gradient stabilized (bottom_solver = 1)
     * CG:      conjugate gradient method (bottom_solver = 2)
     * CABiCG:  (bottom_solver = 3)
     * special: bottom_solver = 4
kraus's avatar
kraus committed
237
     *
frey_m's avatar
frey_m committed
238 239
     * The special bottom solver extends the range of the multigrid coarsening
     * by aggregating coarse grids on the original mesh together and further coarsening.
kraus's avatar
kraus committed
240
     *
frey_m's avatar
frey_m committed
241 242 243 244 245
     * if use_cg == 1 && cg_solver == 0 --> CG
     * if use_cg == 1 && cg_solver == 1 --> BiCG
     * if use_cg == 1 && cg_solver == 2 --> CABiCG
     * if use_cg == 0                   --> CABiCG
     */
kraus's avatar
kraus committed
246 247
    pp_mg.add("bottom_solver", 1);

frey_m's avatar
frey_m committed
248 249 250
    // verbosity of the bottom solver. Higher numbers give more verbosity (doc: cg_verbose)
    amrex::ParmParse pp_cg("cg");
    pp_cg.add("v", 0);
kraus's avatar
kraus committed
251 252


frey_m's avatar
frey_m committed
253
    /* Additional parameters that can't be set by ParmParse:
kraus's avatar
kraus committed
254
     *
frey_m's avatar
frey_m committed
255 256 257 258 259 260 261 262 263
     *      - max_bottom_nlevel = 3 (additional coarsening if you use bottom_solver == 4)
     *      - min_width = 2 (minimum size of grid at coarsest multigrid level)
     *      - max_L0_growth = -1
     */
    amrex::MGT_Solver::def_min_width = 2;
    amrex::MGT_Solver::def_max_L0_growth = -1.0;
}


264 265
double FMGPoissonSolver::solveWithF90_m(const AmrFieldContainer_pt& rho,
                                        const AmrFieldContainer_pt& phi,
frey_m's avatar
frey_m committed
266
                                        const amrex::Vector< AmrFieldContainer_pt > & grad_phi_edge,
267 268 269
                                        const GeomContainer_t& geom,
                                        int baseLevel,
                                        int finestLevel)
frey_m's avatar
frey_m committed
270 271
{
    int nlevs = finestLevel - baseLevel + 1;
kraus's avatar
kraus committed
272

frey_m's avatar
frey_m committed
273 274 275
    GeomContainer_t geom_p(nlevs);
    AmrFieldContainer_pt rho_p(nlevs);
    AmrFieldContainer_pt phi_p(nlevs);
kraus's avatar
kraus committed
276

frey_m's avatar
frey_m committed
277
    for (int ilev = 0; ilev < nlevs; ++ilev) {
frey_m's avatar
frey_m committed
278 279 280
        geom_p[ilev] = geom[ilev + baseLevel];
        rho_p[ilev]  = rho[ilev + baseLevel];
        phi_p[ilev]  = phi[ilev + baseLevel];
frey_m's avatar
frey_m committed
281
    }
kraus's avatar
kraus committed
282

frey_m's avatar
frey_m committed
283 284 285
    //FIXME Refinement ratio
    amrex::IntVect crse_ratio = (baseLevel == 0) ?
        amrex::IntVect::TheZeroVector() : itsAmrObject_mp->refRatio(0);
frey_m's avatar
frey_m committed
286

frey_m's avatar
frey_m committed
287
    amrex::FMultiGrid fmg(geom_p, baseLevel, crse_ratio);
kraus's avatar
kraus committed
288

frey_m's avatar
frey_m committed
289
    if (baseLevel == 0)
290
        fmg.set_bc(bc_m, *phi_p[baseLevel]);
frey_m's avatar
frey_m committed
291
    else
292
        fmg.set_bc(bc_m, *phi_p[baseLevel-1], *phi_p[baseLevel]);
kraus's avatar
kraus committed
293

frey_m's avatar
frey_m committed
294 295
    /* (alpha * a - beta * (del dot b grad)) phi = rho
     * (b = 1)
kraus's avatar
kraus committed
296
     *
frey_m's avatar
frey_m committed
297 298
     * The function call set_const_gravity_coeffs() sets alpha = 0.0
     * and beta = -1 (in MGT_Solver::set_const_gravity_coeffs)
kraus's avatar
kraus committed
299
     *
frey_m's avatar
frey_m committed
300 301 302
     * --> (del dot grad) phi = rho
     */
    fmg.set_const_gravity_coeffs();
kraus's avatar
kraus committed
303

304
    // order of approximation at Dirichlet boundaries
frey_m's avatar
frey_m committed
305 306 307 308
    fmg.set_maxorder(3);

    int always_use_bnorm = 0;
    int need_grad_phi = 1;
309

310
    double residNorm = fmg.solve(phi_p, rho_p, reltol_m, abstol_m, always_use_bnorm, need_grad_phi);
kraus's avatar
kraus committed
311

frey_m's avatar
frey_m committed
312 313 314 315
    for (int ilev = 0; ilev < nlevs; ++ilev) {
        int amr_level = ilev + baseLevel;
        fmg.get_fluxes(grad_phi_edge[amr_level], ilev);
    }
kraus's avatar
kraus committed
316

317
    return residNorm;
frey_m's avatar
frey_m committed
318
}
319 320


frey_m's avatar
frey_m committed
321
/*
322
void FMGPoissonSolver::interpolate_m(AmrScalarFieldContainer_t& phi,
323
                                     const GeomContainer_t& geom,
324 325
                                     double l0norm,
                                     int finestLevel)
326 327 328 329
{
    amrex::PhysBCFunct cphysbc, fphysbc;
    int lo_bc[] = {INT_DIR, INT_DIR, INT_DIR}; // periodic boundaries
    int hi_bc[] = {INT_DIR, INT_DIR, INT_DIR};
frey_m's avatar
frey_m committed
330
    amrex::Vector<amrex::BCRec> bcs(1, amrex::BCRec(lo_bc, hi_bc));
331
    amrex::PCInterp mapper;
kraus's avatar
kraus committed
332

333
    std::unique_ptr<AmrField_t> tmp;
334

335
    for (int lev = 1; lev <= finestLevel; ++lev) {
336 337
        const AmrGrid_t& ba = phi[lev]->boxArray();
        const AmrProcMap_t& dm = phi[lev]->DistributionMap();
338 339
        tmp.reset(new AmrField_t(ba, dm, 1, 0));
        tmp->setVal(0.0);
kraus's avatar
kraus committed
340

341 342
        amrex::InterpFromCoarseLevel(*tmp, 0.0, *phi[lev-1],
                                     0, 0, 1, geom[lev-1], geom[lev],
343
                                     cphysbc, fphysbc,
344
                                     itsAmrObject_mp->refRatio(lev-1),
345
                                     &mapper, bcs);
346 347 348
        phi[lev]->plus(*tmp, 0, 1, 0);
        phi[lev-1]->mult(l0norm, 0, 1);
    }
kraus's avatar
kraus committed
349

350
    phi[finestLevel]->mult(l0norm, 0, 1);
351
}
kraus's avatar
kraus committed
352
*/