AmrMultiGrid.h 24 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5
//
// Class AmrMultiGrid
//   Main class of the AMR Poisson multigrid solver.
//   It implements the multigrid solver described in https://doi.org/10.1016/j.cpc.2019.106912
//
frey_m's avatar
frey_m committed
6
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
// 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/>.
//

23 24 25 26 27 28 29 30
#ifndef AMR_MULTI_GRID_H
#define AMR_MULTI_GRID_H

#include <vector>
#include <memory>

#include <AMReX.H>
#include <AMReX_MultiFab.H>
frey_m's avatar
frey_m committed
31
#include <AMReX_Vector.H>
32 33 34

#include "AmrMultiGridCore.h"

35 36 37
#include "Solvers/AmrPoissonSolver.h"
#include "Amr/AmrBoxLib.h"

38 39
#include "AmrMultiGridLevel.h"

40 41
#include <fstream>

frey_m's avatar
frey_m committed
42 43
#define AMR_MG_TIMER true
#define AMR_MG_WRITE false
44

45
class AmrMultiGrid : public AmrPoissonSolver< AmrBoxLib > {
46 47
    
public:
48 49 50 51 52 53 54 55
    typedef amr::matrix_t         matrix_t;
    typedef amr::vector_t         vector_t;
    typedef amr::multivector_t    mv_t;
    typedef amr::dmap_t           dmap_t;
    typedef amr::comm_t           comm_t;
    typedef amr::local_ordinal_t  lo_t;
    typedef amr::global_ordinal_t go_t;
    typedef amr::scalar_t         scalar_t;
56 57 58
    
    typedef AmrMultiGridLevel<matrix_t, vector_t> AmrMultiGridLevel_t;
    
59 60 61 62 63 64
    typedef AmrMultiGridLevel_t::AmrField_t     AmrField_t;
    typedef AmrMultiGridLevel_t::AmrGeometry_t  AmrGeometry_t;
    typedef AmrMultiGridLevel_t::AmrField_u     AmrField_u;
    typedef AmrMultiGridLevel_t::AmrField_s     AmrField_s;
    typedef AmrMultiGridLevel_t::AmrIntVect_t   AmrIntVect_t;
    typedef AmrMultiGridLevel_t::indices_t      indices_t;
65
    typedef AmrMultiGridLevel_t::coefficients_t coefficients_t;
66 67
    typedef AmrMultiGridLevel_t::umap_t         umap_t;
    typedef AmrMultiGridLevel_t::boundary_t     boundary_t;
68 69 70
    
    typedef BottomSolver<
        Teuchos::RCP<matrix_t>,
frey_m's avatar
frey_m committed
71 72
        Teuchos::RCP<mv_t>,
        AmrMultiGridLevel_t
73 74
    > bsolver_t;
    
frey_m's avatar
frey_m committed
75 76 77 78 79 80 81 82
    typedef BelosBottomSolver<AmrMultiGridLevel_t>      BelosSolver_t;
    typedef Amesos2BottomSolver<AmrMultiGridLevel_t>    Amesos2Solver_t;
    typedef MueLuBottomSolver<AmrMultiGridLevel_t>      MueLuSolver_t;
    
    typedef AmrPreconditioner<matrix_t, AmrMultiGridLevel_t> preconditioner_t;
    
    typedef Ifpack2Preconditioner<AmrMultiGridLevel_t> Ifpack2Preconditioner_t;
    typedef MueLuPreconditioner<AmrMultiGridLevel_t> MueLuPreconditioner_t;
frey_m's avatar
frey_m committed
83
    
84 85 86 87
    typedef amrex::BoxArray boxarray_t;
    typedef amrex::Box box_t;
    typedef amrex::BaseFab<int> basefab_t;
    typedef amrex::FArrayBox farraybox_t;
88
            
89 90
    typedef AmrSmoother::Smoother Smoother;
    
91 92
    typedef amr::Preconditioner Preconditioner;
    
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 123 124 125 126 127 128 129
    /// Supported interpolaters for prolongation operation
    enum Interpolater {
        TRILINEAR = 0,
        LAGRANGE,
        PIECEWISE_CONST
    };
    
    /// Supported bottom solvers
    enum BaseSolver {
        // all Belos
        BICGSTAB,
        MINRES,
        PCPG,
        CG,
        GMRES,
        STOCHASTIC_CG,
        RECYCLING_CG,
        RECYCLING_GMRES
        // all Amesos2
#ifdef HAVE_AMESOS2_KLU2
        , KLU2
#endif
#ifdef HAVE_AMESOS2_SUPERLU
        , SUPERLU
#endif
#ifdef HAVE_AMESOS2_UMFPACK  
        , UMFPACK
#endif
#ifdef HAVE_AMESOS2_PARDISO_MKL
        , PARDISO_MKL
#endif
#ifdef HAVE_AMESOS2_MUMPS
        , MUMPS
#endif
#ifdef HAVE_AMESOS2_LAPACK
        , LAPACK
#endif
frey_m's avatar
frey_m committed
130 131
        // all MueLu
        , SA
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
        // add others ...
    };
    
    /// Supported physical boundaries
    enum Boundary {
        DIRICHLET = 0,
        OPEN,
        PERIODIC
    };
    
    /// Supported convergence criteria
    enum Norm {
        L1,
        L2,
        LINF
    };
    
public:
    
    /*!
     * Instantiation used in Structure/FieldSolver.cpp
153
     * @param itsAmrObject_p has information about refinemen ratios, etc.
154 155
     * @param bsolver bottom solver
     * @param prec preconditioner for bottom solver
156
     * @param rebalance of preconditioner (SA only)
157 158 159 160 161 162 163 164
     * @param bcx boundary condition in x
     * @param bcy boundary condition in y
     * @param bcz boundary condition in z
     * @param smoother for level solution
     * @param nSweeps when smoothing
     * @param interp interpolater between levels
     * @param norm for convergence criteria
     */
165 166
    AmrMultiGrid(AmrBoxLib* itsAmrObject_p,
                 const std::string& bsolver,
167
                 const std::string& prec,
168
                 const bool& rebalance,
frey_m's avatar
frey_m committed
169
                 const std::string& reuse,
170 171 172 173 174 175 176 177
                 const std::string& bcx,
                 const std::string& bcy,
                 const std::string& bcz,
                 const std::string& smoother,
                 const std::size_t& nSweeps,
                 const std::string& interp,
                 const std::string& norm);
    
178 179 180 181 182 183 184 185 186 187
    /*!
     * Used in OPAL
     * 
     * @param rho right-hand side charge density on grid [C / m]
     * @param phi electrostatic potential (unknown) [V]
     * @param efield electric field [V / m]
     * @param baseLevel for solve
     * @param finestLevel for solve
     * @param prevAsGuess use of previous solution as initial guess
     */
188 189 190
    void solve(AmrScalarFieldContainer_t &rho,
               AmrScalarFieldContainer_t &phi,
               AmrVectorFieldContainer_t &efield,
frey_m's avatar
frey_m committed
191 192 193
               unsigned short baseLevel,
               unsigned short finestLevel,
               bool prevAsGuess = true);
194
    
195 196 197 198
    /*!
     * Specify the number of smoothing steps
     * @param nSweeps for each smoothing step
     */
199
    void setNumberOfSweeps(const std::size_t& nSweeps);
200
    
201 202 203 204 205
    /*!
     * Specify the maximum number of iterations
     * @param maxiter \f$ [0, \infty[ \f$
     */
    void setMaxNumberOfIterations(const std::size_t& maxiter);
206 207 208 209 210
    
    /*!
     * Obtain some convergence info
     * @returns the number of iterations till convergence
     */
211
    std::size_t getNumIters();
212
    
213 214 215 216 217
    /*!
     * Obtain the residual norm of a level
     * @param level for which error is requested
     * @returns the norm of the residual
     */
218
    scalar_t getLevelResidualNorm(lo_t level);
219
    
220 221 222 223
    /*!
     * Enable solver info dumping into SDDS file
     */
    void setVerbose(bool verbose);
224
    
225 226 227 228 229 230 231 232 233 234 235 236
    double getXRangeMin(unsigned short level = 0);
    double getXRangeMax(unsigned short level = 0);
    double getYRangeMin(unsigned short level = 0);
    double getYRangeMax(unsigned short level = 0);
    double getZRangeMin(unsigned short level = 0);
    double getZRangeMax(unsigned short level = 0);
    
    /**
     * Print information abour tolerances.
     * @param os output stream where to write to
     */
    Inform &print(Inform &os) const;
237 238 239 240 241
    
private:
    
    /*!
     * Instantiate boundary object
242 243
     * @param bc boundary conditions
     * @precondition length must be equal to AMREX_SPACEDIM
244
     */
245
    void initPhysicalBoundary_m(const Boundary* bc);
246 247 248 249 250
    
    /*!
     * Instantiate all levels and set boundary conditions
     * @param rho is the charge density
     * @param geom is the geometry
frey_m's avatar
frey_m committed
251
     * @param regrid was performed
252
     */
frey_m's avatar
frey_m committed
253 254
    void initLevels_m(const amrex::Vector<AmrField_u>& rho,
                      const amrex::Vector<AmrGeometry_t>& geom,
frey_m's avatar
frey_m committed
255
                      bool regrid);
256
    
257 258 259 260 261
    /*!
     * Clear masks (required to build matrices) no longer needed.
     */
    void clearMasks_m();
    
262 263
    /*!
     * Reset potential to zero (currently)
frey_m's avatar
frey_m committed
264
     * @param reset solution to initial guess (zero)
265
     */
frey_m's avatar
frey_m committed
266
    void initGuess_m(bool reset);
267 268 269
    
    /*!
     * Actual solve.
270
     * @returns the the max. residual
271
     */
272
    scalar_t iterate_m();
frey_m's avatar
frey_m committed
273
    
frey_m's avatar
frey_m committed
274 275 276 277
    /*!
     * Compute norms / level and check convergence
     * @returns true if converged
     */
frey_m's avatar
frey_m committed
278 279
    bool isConverged_m(std::vector<scalar_t>& rhsNorms,
                       std::vector<scalar_t>& resNorms);
frey_m's avatar
frey_m committed
280
    
281 282 283 284 285 286 287
    /*!
     * Compute composite residual of a level
     * @param r is the residual to compute
     * @param b is the right-hand side
     * @param x is the left-hand side
     * @param level to solve for
     */
288 289
    void residual_m(const lo_t& level,
                    Teuchos::RCP<vector_t>& r,
290
                    const Teuchos::RCP<vector_t>& b,
291
                    const Teuchos::RCP<vector_t>& x);
292 293 294 295 296
    
    /*!
     * Recursive call.
     * @param level to relax
     */
297
    void relax_m(const lo_t& level);
298 299 300 301 302 303 304 305 306
    
    /*!
     * Compute the residual of a level without considering refined level.
     * @param result is computed
     * @param rhs is the right-hand side
     * @param crs_rhs is the coarse right-hand side for internal boundary
     * @param b is the left-hand side
     * @param level to solver for
     */
307 308
    void residual_no_fine_m(const lo_t& level,
                            Teuchos::RCP<vector_t>& result,
309 310
                            const Teuchos::RCP<vector_t>& rhs,
                            const Teuchos::RCP<vector_t>& crs_rhs,
311
                            const Teuchos::RCP<vector_t>& b);
frey_m's avatar
frey_m committed
312
    
frey_m's avatar
frey_m committed
313
#if AMR_MG_WRITE
314
    /*!
frey_m's avatar
frey_m committed
315
     * Dumps the residual norm per level into a file (for each iteration).
316
     */
frey_m's avatar
frey_m committed
317 318
    void writeResidualNorm_m();
#endif
319 320 321 322 323 324
    
    /*!
     * Vector norm computation.
     * @param x is the vector for which we compute the norm
     * @returns the evaluated norm of a level
     */
325
    scalar_t evalNorm_m(const Teuchos::RCP<const vector_t>& x);
326 327 328
    
    /*!
     * Initial convergence criteria values.
frey_m's avatar
frey_m committed
329 330
     * @param rhsNorms per level of right-hand side (is filled)
     * @param resNorms per level of residual (is filled)
331
     */
frey_m's avatar
frey_m committed
332 333
    void initResidual_m(std::vector<scalar_t>& rhsNorms,
                        std::vector<scalar_t>& resNorms);
334 335 336 337
    
    /*!
     * @param efield to compute
     */
338
    void computeEfield_m(AmrVectorFieldContainer_t& efield);
339 340 341 342 343
    
    /*!
     * Build all matrices and vectors, i.e. AMReX to Trilinos
     * @param matrices if we need to build matrices as well or only vectors.
     */
frey_m's avatar
frey_m committed
344 345
    void setup_m(const amrex::Vector<AmrField_u>& rho,
                 const amrex::Vector<AmrField_u>& phi,
346
                 const bool& matrices = true);
347
    
348 349 350
    /*!
     * Build all matrices and vectors needed for single-level computation
     */
frey_m's avatar
frey_m committed
351 352
    void buildSingleLevel_m(const amrex::Vector<AmrField_u>& rho,
                            const amrex::Vector<AmrField_u>& phi,
353 354 355 356 357
                            const bool& matrices = true);
    
    /*!
     * Build all matrices and vectors needed for multi-level computation
     */
frey_m's avatar
frey_m committed
358 359
    void buildMultiLevel_m(const amrex::Vector<AmrField_u>& rho,
                           const amrex::Vector<AmrField_u>& phi,
360 361
                           const bool& matrices = true);

362 363 364 365 366
    /*!
     * Set matrix and vector pointer
     * @param level for which we fill matrix + vector
     * @param matrices if we need to set matrices as well or only vectors.
     */
367
    void open_m(const lo_t& level, const bool& matrices);
368 369 370 371 372 373
    
    /*!
     * Call fill complete
     * @param level for which we filled matrix
     * @param matrices if we set matrices as well.
     */
374
    void close_m(const lo_t& level, const bool& matrices);
375 376 377 378 379 380 381 382 383 384 385
    
    /*!
     * Build the Poisson matrix for a level assuming no finer level (i.e. the whole fine mesh
     * is taken into account).
     * It takes care of physical boundaries (i.e. mesh boundary).
     * Internal boundaries (i.e. boundaries due to crse-fine interfaces) are treated by the
     * boundary matrix.
     * @param iv is the current cell
     * @param mfab is the mask (internal cell, boundary cell, ...) of that level
     * @param level for which we build the Poisson matrix
     */
386 387 388 389 390
    void buildNoFinePoissonMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& mfab,
                                    const scalar_t* invdx2);
391 392 393 394 395 396 397 398 399
    
    /*!
     * Build the Poisson matrix for a level that got refined (it does not take the covered
     * cells (covered by fine cells) into account). The finest level does not build such a matrix.
     * It takes care of physical boundaries (i.e. mesh boundary).
     * Internal boundaries (i.e. boundaries due to crse-fine interfaces) are treated by the
     * boundary matrix.
     * @param iv is the current cell
     * @param mfab is the mask (internal cell, boundary cell, ...) of that level
400
     * @param rfab is the mask between levels
401 402
     * @param level for which we build the special Poisson matrix
     */
403 404 405
    void buildCompositePoissonMatrix_m(const lo_t& level,
                                       const go_t& gidx,
                                       const AmrIntVect_t& iv,
406
                                       const basefab_t& mfab,
407 408 409
                                       const basefab_t& rfab,
                                       const basefab_t& cfab,
                                       const scalar_t* invdx2);
410 411 412 413 414 415 416 417
    
    /*!
     * Build a matrix that averages down the data of the fine cells down to the
     * corresponding coarse cells. The base level does not build such a matrix.
     * \f[
     *      x^{(l)} = R\cdot x^{(l+1)}
     * \f]
     * @param iv is the current cell
418
     * @param rfab is the mask between levels
419 420
     * @param level for which to build restriction matrix
     */
421 422 423 424 425 426 427
    void buildRestrictionMatrix_m(const lo_t& level,
                                  const go_t& gidx,
                                  const AmrIntVect_t& iv,
                                  D_DECL(const go_t& ii,
                                         const go_t& jj,
                                         const go_t& kk),
                                  const basefab_t& rfab);
428 429 430 431 432 433 434 435 436 437 438 439
    
    /*!
     * Interpolate data from coarse cells to appropriate refined cells. The interpolation
     * scheme is allowed only to have local cells in the stencil, i.e.
     * 2D --> 4, 3D --> 8.
     * \f[
     *      x^{(l)} = I\cdot x^{(l-1)}
     * \f]
     * @param iv is the current cell
     * @param level for which to build the interpolation matrix. The finest level
     * does not build such a matrix.
     */
440 441 442 443
    void buildInterpolationMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& cfab);
444 445 446 447 448 449 450 451 452
    
    /*!
     * The boundary values at the crse-fine-interface need to be taken into account properly.
     * This matrix is used to update the fine boundary values from the coarse values, i.e.
     * \f[
     *      x^{(l)} = B_{crse}\cdot x^{(l-1)}
     * \f]
     * Dirichlet boundary condition
     * @param iv is the current cell
453
     * @param mfab is the mask (internal cell, boundary cell, ...) of that level
454 455 456
     * @param cells all fine cells that are at the crse-fine interface
     * @param level the base level is omitted
     */
457 458 459 460 461 462
    void buildCrseBoundaryMatrix_m(const lo_t& level,
                                   const go_t& gidx,
                                   const AmrIntVect_t& iv,
                                   const basefab_t& mfab,
                                   const basefab_t& cfab,
                                   const scalar_t* invdx2);
463 464 465 466 467 468 469 470 471 472 473 474 475 476
    
    /*!
     * The boundary values at the crse-fine-interface need to be taken into account properly.
     * This matrix is used to update the coarse boundary values from fine values, i.e.
     * \f[
     *      x^{(l)} = B_{fine}\cdot x^{(l+1)}
     * \f]
     * Dirichlet boundary condition. Flux matching.
     * @param iv is the current cell
     * @param cells all coarse cells that are at the crse-fine interface but are
     * not refined
     * @param crse_fine_ba coarse cells that got refined
     * @param level the finest level is omitted
     */
477 478 479 480 481 482
    void buildFineBoundaryMatrix_m(const lo_t& level,
                                   const go_t& gidx,
                                   const AmrIntVect_t& iv,
                                   const basefab_t& mfab,
                                   const basefab_t& rfab,
                                   const basefab_t& cfab);
483 484 485 486 487 488
    
    /*!
     * Copy data from AMReX to Trilinos
     * @param rho is the charge density
     * @param level for which to copy
     */
489 490
    void buildDensityVector_m(const lo_t& level,
                              const AmrField_t& rho);
491 492 493 494 495 496
    
    /*!
     * Copy data from AMReX to Trilinos
     * @param phi is the potential
     * @param level for which to copy
     */
497 498
    void buildPotentialVector_m(const lo_t& level,
                                const AmrField_t& phi);
499 500 501 502 503 504 505
    
    /*!
     * Gradient matrix is used to compute the electric field
     * @param iv is the current cell
     * @param mfab is the mask (internal cell, boundary cell, ...) of that level
     * @param level for which to compute
     */
506 507 508 509 510
    void buildGradientMatrix_m(const lo_t& level,
                               const go_t& gidx,
                               const AmrIntVect_t& iv,
                               const basefab_t& mfab,
                               const scalar_t* invdx);
511 512 513 514 515 516 517 518
    
    /*!
     * Data transfer from AMReX to Trilinos.
     * @param mf is the multifab of a level
     * @param comp component to copy
     * @param mv is the vector to be filled
     * @param level where to perform
     */
519 520 521 522
    void amrex2trilinos_m(const lo_t& level,
                          const lo_t& comp,
                          const AmrField_t& mf,
                          Teuchos::RCP<vector_t>& mv);
523 524 525
    
    /*!
     * Data transfer from Trilinos to AMReX.
526
     * @param level to copy
527
     * @param comp component to copy
528
     * @param mf is the multifab to be filled
529 530
     * @param mv is the corresponding Trilinos vector
     */
531 532
    void trilinos2amrex_m(const lo_t& level,
                          const lo_t& comp,
533 534
                          AmrField_t& mf,
                          const Teuchos::RCP<vector_t>& mv);
535 536
    
    /*!
537 538 539
     * Some indices might occur several times due to boundary conditions, etc. We
     * avoid this by filling a map and then copying the data to a vector for filling
     * the matrices. The map gets cleared inside the function.
540 541 542
     * @param indices in matrix
     * @param values are the coefficients
     */
543 544
    void map2vector_m(umap_t& map, indices_t& indices,
                      coefficients_t& values);
545 546
    
    /*!
547
     * Perform one smoothing step
548 549 550 551
     * @param e error to update (left-hand side)
     * @param r residual (right-hand side)
     * @param level on which to relax
     */
552 553 554
    void smooth_m(const lo_t& level,
                  Teuchos::RCP<vector_t>& e,
                  Teuchos::RCP<vector_t>& r);
555 556 557 558 559
    
    /*!
     * Restrict coarse level residual based on fine level residual
     * @param level to restrict
     */
560
    void restrict_m(const lo_t& level);
561 562 563 564 565
    
    /*!
     * Update error of fine level based on error of coarse level
     * @param level to update
     */
566
    void prolongate_m(const lo_t& level);
567 568 569 570 571
    
    /*!
     * Average data from fine level to coarse
     * @param level finest level is omitted
     */
572
    void averageDown_m(const lo_t& level);
573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588
    
    /*!
     * Instantiate interpolater
     * @param interp interpolater type
     */
    void initInterpolater_m(const Interpolater& interp);
    
    /*!
     * Instantiate interface interpolater
     * @param interface handler
     */
    void initCrseFineInterp_m(const Interpolater& interface);
    
    /*!
     * Instantiate a bottom solver
     * @param solver type
frey_m's avatar
frey_m committed
589 590
     * @param rebalance solver (SA only)
     * @param reuse types of SA hierarchy
591
     */
592
    void initBaseSolver_m(const BaseSolver& solver,
frey_m's avatar
frey_m committed
593 594
                          const bool& rebalance,
                          const std::string& reuse);
595 596 597 598
    
    /*!
     * Instantiate a preconditioner for the bottom solver
     * @param precond type
599
     * @param rebalance preconditioner (SA only)
frey_m's avatar
frey_m committed
600
     * @param reuse types of SA hierarchy
601
     */
602
    void initPrec_m(const Preconditioner& prec,
frey_m's avatar
frey_m committed
603 604
                    const bool& rebalance,
                    const std::string& reuse);
605
    
606 607 608 609 610 611 612 613 614 615 616 617
    /*!
     * Convertstring to enum Boundary
     * @param bc boundary condition
     */
    Boundary convertToEnumBoundary_m(const std::string& bc);
    
    /*!
     * Converts string to enum Interpolater
     * @param interp interpolater
     */
    Interpolater convertToEnumInterpolater_m(const std::string& interp);
    
618 619 620 621
    /*!
     * Converts string to enum BaseSolver
     * @param bsolver bottom solver
     */
622
    BaseSolver convertToEnumBaseSolver_m(const std::string& bsolver);
623 624 625 626 627
    
    /*!
     * Converts string to enum Preconditioner
     * @param prec preconditioner
     */
628 629 630 631 632 633 634
    Preconditioner convertToEnumPreconditioner_m(const std::string& prec);
    
    /*!
     * Converts string to enum Smoother
     * @param smoother of level solution
     */
    Smoother convertToEnumSmoother_m(const std::string& smoother);
635
    
636 637 638 639 640
    /*!
     * Converts string to enum Norm
     * @param norm either L1, L2, LInf
     */
    Norm convertToEnumNorm_m(const std::string& norm);
641
    
642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
    /*!
     * SDDS header is written by root core
     * @param outfile output stream
     */
    void writeSDDSHeader_m(std::ofstream& outfile);
    
    /*!
     * SDDS data write (done by root core)
     * @param error to write
     */
    void writeSDDSData_m(const scalar_t& error);
    
#if AMR_MG_TIMER
    /*!
     * Create timers
     */
    void initTimer_m();
#endif
    
661 662 663 664 665 666 667 668 669 670 671
private:
    Teuchos::RCP<comm_t> comm_mp;       ///< communicator
    Teuchos::RCP<amr::node_t> node_mp;  ///< kokkos node
    
    /// interpolater without coarse-fine interface
    std::unique_ptr<AmrInterpolater<AmrMultiGridLevel_t> > interp_mp;
    
    /// interpolater for coarse-fine interface
    std::unique_ptr<AmrInterpolater<AmrMultiGridLevel_t> > interface_mp;
    
    std::size_t nIter_m;            ///< number of iterations till convergence
672
    std::size_t bIter_m;            ///< number of iterations of bottom solver
673
    std::size_t maxiter_m;          ///< maximum number of iterations allowed
674 675 676 677 678 679 680
    std::size_t nSweeps_m;          ///< number of smoothing iterations
    Smoother smootherType_m;        ///< type of smoother
    
    /// container for levels
    std::vector<std::unique_ptr<AmrMultiGridLevel_t > > mglevel_m;
    
    /// bottom solver
681
    std::shared_ptr<bsolver_t> solver_mp;
682 683 684 685
    
    /// error smoother
    std::vector<std::shared_ptr<AmrSmoother> > smoother_m;
    
686 687 688
    /// preconditioner for bottom solver
    std::shared_ptr<preconditioner_t> prec_mp;
    
689 690 691 692
    int lbase_m;            ///< base level (currently only 0 supported)
    int lfine_m;            ///< fineste level
    int nlevel_m;           ///< number of levelss
    
693 694
    boundary_t bc_m[AMREX_SPACEDIM];    ///< boundary conditions
    int nBcPoints_m;                    ///< maximum number of stencils points for BC
695 696
    
    Norm norm_m;            ///< norm for convergence criteria (l1, l2, linf)
697
    std::string snorm_m;    ///< norm for convergence criteria
698
    
699
    const scalar_t eps_m;   ///< rhs scale for convergence
700
    
701 702 703 704
    bool verbose_m;                 ///< If true, a SDDS file is written
    std::string fname_m;            ///< SDDS filename
    std::ios_base::openmode flag_m; ///< std::ios::out or std::ios::app
    
705 706 707 708 709
#if AMR_MG_TIMER
    IpplTimings::TimerRef buildTimer_m;         ///< timer for matrix and vector construction
    IpplTimings::TimerRef restrictTimer_m;      ///< timer for restriction operation
    IpplTimings::TimerRef smoothTimer_m;        ///< timer for all smoothing steps
    IpplTimings::TimerRef interpTimer_m;        ///< prolongation timer
710 711
    IpplTimings::TimerRef efieldTimer_m;        ///< timer for e-field calculation + copy
    IpplTimings::TimerRef averageTimer_m;       ///> timer for average down process
712
    IpplTimings::TimerRef bottomTimer_m;        ///< bottom solver timer
713
    IpplTimings::TimerRef dumpTimer_m;          ///< write SDDS file timer
714 715 716
#endif
};

717 718 719 720 721

inline Inform &operator<<(Inform &os, const AmrMultiGrid &fs) {
    return fs.print(os);
}

722
#endif