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
// 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 27 28 29
#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
30
#include <AMReX_Vector.H>
31 32 33

#include "AmrMultiGridCore.h"

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

37 38
#include "AmrMultiGridLevel.h"

39 40
#include <fstream>

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

44
class AmrMultiGrid : public AmrPoissonSolver< AmrBoxLib > {
45 46
    
public:
47 48 49 50 51 52 53 54
    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;
55 56 57
    
    typedef AmrMultiGridLevel<matrix_t, vector_t> AmrMultiGridLevel_t;
    
58 59 60 61 62 63
    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;
64
    typedef AmrMultiGridLevel_t::coefficients_t coefficients_t;
65 66
    typedef AmrMultiGridLevel_t::umap_t         umap_t;
    typedef AmrMultiGridLevel_t::boundary_t     boundary_t;
67 68 69
    
    typedef BottomSolver<
        Teuchos::RCP<matrix_t>,
frey_m's avatar
frey_m committed
70 71
        Teuchos::RCP<mv_t>,
        AmrMultiGridLevel_t
72 73
    > bsolver_t;
    
frey_m's avatar
frey_m committed
74 75 76 77 78 79 80 81
    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
82
    
83 84 85 86
    typedef amrex::BoxArray boxarray_t;
    typedef amrex::Box box_t;
    typedef amrex::BaseFab<int> basefab_t;
    typedef amrex::FArrayBox farraybox_t;
87
            
88 89
    typedef AmrSmoother::Smoother Smoother;
    
90 91
    typedef amr::Preconditioner Preconditioner;
    
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 123 124 125 126 127 128
    /// 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
129 130
        // all MueLu
        , SA
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
        // 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
152
     * @param itsAmrObject_p has information about refinemen ratios, etc.
153 154
     * @param bsolver bottom solver
     * @param prec preconditioner for bottom solver
155
     * @param rebalance of preconditioner (SA only)
156 157 158 159 160 161 162 163
     * @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
     */
164 165
    AmrMultiGrid(AmrBoxLib* itsAmrObject_p,
                 const std::string& bsolver,
166
                 const std::string& prec,
167
                 const bool& rebalance,
frey_m's avatar
frey_m committed
168
                 const std::string& reuse,
169 170 171 172 173 174 175 176
                 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);
    
177 178 179 180 181 182 183 184 185 186
    /*!
     * 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
     */
187 188 189
    void solve(AmrScalarFieldContainer_t &rho,
               AmrScalarFieldContainer_t &phi,
               AmrVectorFieldContainer_t &efield,
frey_m's avatar
frey_m committed
190 191 192
               unsigned short baseLevel,
               unsigned short finestLevel,
               bool prevAsGuess = true);
193
    
194 195 196 197
    /*!
     * Specify the number of smoothing steps
     * @param nSweeps for each smoothing step
     */
198
    void setNumberOfSweeps(const std::size_t& nSweeps);
199
    
200 201 202 203 204
    /*!
     * Specify the maximum number of iterations
     * @param maxiter \f$ [0, \infty[ \f$
     */
    void setMaxNumberOfIterations(const std::size_t& maxiter);
205 206 207 208 209
    
    /*!
     * Obtain some convergence info
     * @returns the number of iterations till convergence
     */
210
    std::size_t getNumIters();
211
    
212 213 214 215 216
    /*!
     * Obtain the residual norm of a level
     * @param level for which error is requested
     * @returns the norm of the residual
     */
217
    scalar_t getLevelResidualNorm(lo_t level);
218
    
219 220 221 222
    /*!
     * Enable solver info dumping into SDDS file
     */
    void setVerbose(bool verbose);
223
    
224 225 226 227 228 229 230 231 232 233 234 235
    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;
236 237 238 239 240
    
private:
    
    /*!
     * Instantiate boundary object
241 242
     * @param bc boundary conditions
     * @precondition length must be equal to AMREX_SPACEDIM
243
     */
244
    void initPhysicalBoundary_m(const Boundary* bc);
245 246 247 248 249
    
    /*!
     * 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
250
     * @param regrid was performed
251
     */
frey_m's avatar
frey_m committed
252 253
    void initLevels_m(const amrex::Vector<AmrField_u>& rho,
                      const amrex::Vector<AmrGeometry_t>& geom,
frey_m's avatar
frey_m committed
254
                      bool regrid);
255
    
256 257 258 259 260
    /*!
     * Clear masks (required to build matrices) no longer needed.
     */
    void clearMasks_m();
    
261 262
    /*!
     * Reset potential to zero (currently)
frey_m's avatar
frey_m committed
263
     * @param reset solution to initial guess (zero)
264
     */
frey_m's avatar
frey_m committed
265
    void initGuess_m(bool reset);
266 267 268
    
    /*!
     * Actual solve.
269
     * @returns the the max. residual
270
     */
271
    scalar_t iterate_m();
frey_m's avatar
frey_m committed
272
    
frey_m's avatar
frey_m committed
273 274 275 276
    /*!
     * Compute norms / level and check convergence
     * @returns true if converged
     */
frey_m's avatar
frey_m committed
277 278
    bool isConverged_m(std::vector<scalar_t>& rhsNorms,
                       std::vector<scalar_t>& resNorms);
frey_m's avatar
frey_m committed
279
    
280 281 282 283 284 285 286
    /*!
     * 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
     */
287 288
    void residual_m(const lo_t& level,
                    Teuchos::RCP<vector_t>& r,
289
                    const Teuchos::RCP<vector_t>& b,
290
                    const Teuchos::RCP<vector_t>& x);
291 292 293 294 295
    
    /*!
     * Recursive call.
     * @param level to relax
     */
296
    void relax_m(const lo_t& level);
297 298 299 300 301 302 303 304 305
    
    /*!
     * 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
     */
306 307
    void residual_no_fine_m(const lo_t& level,
                            Teuchos::RCP<vector_t>& result,
308 309
                            const Teuchos::RCP<vector_t>& rhs,
                            const Teuchos::RCP<vector_t>& crs_rhs,
310
                            const Teuchos::RCP<vector_t>& b);
frey_m's avatar
frey_m committed
311
    
frey_m's avatar
frey_m committed
312
#if AMR_MG_WRITE
313
    /*!
frey_m's avatar
frey_m committed
314
     * Dumps the residual norm per level into a file (for each iteration).
315
     */
frey_m's avatar
frey_m committed
316 317
    void writeResidualNorm_m();
#endif
318 319 320 321 322 323
    
    /*!
     * Vector norm computation.
     * @param x is the vector for which we compute the norm
     * @returns the evaluated norm of a level
     */
324
    scalar_t evalNorm_m(const Teuchos::RCP<const vector_t>& x);
325 326 327
    
    /*!
     * Initial convergence criteria values.
frey_m's avatar
frey_m committed
328 329
     * @param rhsNorms per level of right-hand side (is filled)
     * @param resNorms per level of residual (is filled)
330
     */
frey_m's avatar
frey_m committed
331 332
    void initResidual_m(std::vector<scalar_t>& rhsNorms,
                        std::vector<scalar_t>& resNorms);
333 334 335 336
    
    /*!
     * @param efield to compute
     */
337
    void computeEfield_m(AmrVectorFieldContainer_t& efield);
338 339 340 341 342
    
    /*!
     * 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
343 344
    void setup_m(const amrex::Vector<AmrField_u>& rho,
                 const amrex::Vector<AmrField_u>& phi,
345
                 const bool& matrices = true);
346
    
347 348 349
    /*!
     * Build all matrices and vectors needed for single-level computation
     */
frey_m's avatar
frey_m committed
350 351
    void buildSingleLevel_m(const amrex::Vector<AmrField_u>& rho,
                            const amrex::Vector<AmrField_u>& phi,
352 353 354 355 356
                            const bool& matrices = true);
    
    /*!
     * Build all matrices and vectors needed for multi-level computation
     */
frey_m's avatar
frey_m committed
357 358
    void buildMultiLevel_m(const amrex::Vector<AmrField_u>& rho,
                           const amrex::Vector<AmrField_u>& phi,
359 360
                           const bool& matrices = true);

361 362 363 364 365
    /*!
     * 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.
     */
366
    void open_m(const lo_t& level, const bool& matrices);
367 368 369 370 371 372
    
    /*!
     * Call fill complete
     * @param level for which we filled matrix
     * @param matrices if we set matrices as well.
     */
373
    void close_m(const lo_t& level, const bool& matrices);
374 375 376 377 378 379 380 381 382 383 384
    
    /*!
     * 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
     */
385 386 387 388 389
    void buildNoFinePoissonMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& mfab,
                                    const scalar_t* invdx2);
390 391 392 393 394 395 396 397 398
    
    /*!
     * 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
399
     * @param rfab is the mask between levels
400 401
     * @param level for which we build the special Poisson matrix
     */
402 403 404
    void buildCompositePoissonMatrix_m(const lo_t& level,
                                       const go_t& gidx,
                                       const AmrIntVect_t& iv,
405
                                       const basefab_t& mfab,
406 407 408
                                       const basefab_t& rfab,
                                       const basefab_t& cfab,
                                       const scalar_t* invdx2);
409 410 411 412 413 414 415 416
    
    /*!
     * 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
417
     * @param rfab is the mask between levels
418 419
     * @param level for which to build restriction matrix
     */
420 421 422 423 424 425 426
    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);
427 428 429 430 431 432 433 434 435 436 437 438
    
    /*!
     * 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.
     */
439 440 441 442
    void buildInterpolationMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& cfab);
443 444 445 446 447 448 449 450 451
    
    /*!
     * 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
452
     * @param mfab is the mask (internal cell, boundary cell, ...) of that level
453 454 455
     * @param cells all fine cells that are at the crse-fine interface
     * @param level the base level is omitted
     */
456 457 458 459 460 461
    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);
462 463 464 465 466 467 468 469 470 471 472 473 474 475
    
    /*!
     * 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
     */
476 477 478 479 480 481
    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);
482 483 484 485 486 487
    
    /*!
     * Copy data from AMReX to Trilinos
     * @param rho is the charge density
     * @param level for which to copy
     */
488 489
    void buildDensityVector_m(const lo_t& level,
                              const AmrField_t& rho);
490 491 492 493 494 495
    
    /*!
     * Copy data from AMReX to Trilinos
     * @param phi is the potential
     * @param level for which to copy
     */
496 497
    void buildPotentialVector_m(const lo_t& level,
                                const AmrField_t& phi);
498 499 500 501 502 503 504
    
    /*!
     * 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
     */
505 506 507 508 509
    void buildGradientMatrix_m(const lo_t& level,
                               const go_t& gidx,
                               const AmrIntVect_t& iv,
                               const basefab_t& mfab,
                               const scalar_t* invdx);
510 511 512 513 514 515 516 517
    
    /*!
     * 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
     */
518 519 520 521
    void amrex2trilinos_m(const lo_t& level,
                          const lo_t& comp,
                          const AmrField_t& mf,
                          Teuchos::RCP<vector_t>& mv);
522 523 524
    
    /*!
     * Data transfer from Trilinos to AMReX.
525
     * @param level to copy
526
     * @param comp component to copy
527
     * @param mf is the multifab to be filled
528 529
     * @param mv is the corresponding Trilinos vector
     */
530 531
    void trilinos2amrex_m(const lo_t& level,
                          const lo_t& comp,
532 533
                          AmrField_t& mf,
                          const Teuchos::RCP<vector_t>& mv);
534 535
    
    /*!
536 537 538
     * 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.
539 540 541
     * @param indices in matrix
     * @param values are the coefficients
     */
542 543
    void map2vector_m(umap_t& map, indices_t& indices,
                      coefficients_t& values);
544 545
    
    /*!
546
     * Perform one smoothing step
547 548 549 550
     * @param e error to update (left-hand side)
     * @param r residual (right-hand side)
     * @param level on which to relax
     */
551 552 553
    void smooth_m(const lo_t& level,
                  Teuchos::RCP<vector_t>& e,
                  Teuchos::RCP<vector_t>& r);
554 555 556 557 558
    
    /*!
     * Restrict coarse level residual based on fine level residual
     * @param level to restrict
     */
559
    void restrict_m(const lo_t& level);
560 561 562 563 564
    
    /*!
     * Update error of fine level based on error of coarse level
     * @param level to update
     */
565
    void prolongate_m(const lo_t& level);
566 567 568 569 570
    
    /*!
     * Average data from fine level to coarse
     * @param level finest level is omitted
     */
571
    void averageDown_m(const lo_t& level);
572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
    
    /*!
     * 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
588 589
     * @param rebalance solver (SA only)
     * @param reuse types of SA hierarchy
590
     */
591
    void initBaseSolver_m(const BaseSolver& solver,
frey_m's avatar
frey_m committed
592 593
                          const bool& rebalance,
                          const std::string& reuse);
594 595 596 597
    
    /*!
     * Instantiate a preconditioner for the bottom solver
     * @param precond type
598
     * @param rebalance preconditioner (SA only)
frey_m's avatar
frey_m committed
599
     * @param reuse types of SA hierarchy
600
     */
601
    void initPrec_m(const Preconditioner& prec,
frey_m's avatar
frey_m committed
602 603
                    const bool& rebalance,
                    const std::string& reuse);
604
    
605 606 607 608 609 610 611 612 613 614 615 616
    /*!
     * 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);
    
617 618 619 620
    /*!
     * Converts string to enum BaseSolver
     * @param bsolver bottom solver
     */
621
    BaseSolver convertToEnumBaseSolver_m(const std::string& bsolver);
622 623 624 625 626
    
    /*!
     * Converts string to enum Preconditioner
     * @param prec preconditioner
     */
627 628 629 630 631 632 633
    Preconditioner convertToEnumPreconditioner_m(const std::string& prec);
    
    /*!
     * Converts string to enum Smoother
     * @param smoother of level solution
     */
    Smoother convertToEnumSmoother_m(const std::string& smoother);
634
    
635 636 637 638 639
    /*!
     * Converts string to enum Norm
     * @param norm either L1, L2, LInf
     */
    Norm convertToEnumNorm_m(const std::string& norm);
640
    
641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659
    /*!
     * 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
    
660 661 662 663 664 665 666 667 668 669 670
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
671
    std::size_t bIter_m;            ///< number of iterations of bottom solver
672
    std::size_t maxiter_m;          ///< maximum number of iterations allowed
673 674 675 676 677 678 679
    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
680
    std::shared_ptr<bsolver_t> solver_mp;
681 682 683 684
    
    /// error smoother
    std::vector<std::shared_ptr<AmrSmoother> > smoother_m;
    
685 686 687
    /// preconditioner for bottom solver
    std::shared_ptr<preconditioner_t> prec_mp;
    
688 689 690 691
    int lbase_m;            ///< base level (currently only 0 supported)
    int lfine_m;            ///< fineste level
    int nlevel_m;           ///< number of levelss
    
692 693
    boundary_t bc_m[AMREX_SPACEDIM];    ///< boundary conditions
    int nBcPoints_m;                    ///< maximum number of stencils points for BC
694 695
    
    Norm norm_m;            ///< norm for convergence criteria (l1, l2, linf)
696
    std::string snorm_m;    ///< norm for convergence criteria
697
    
698
    const scalar_t eps_m;   ///< rhs scale for convergence
699
    
700 701 702 703
    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
    
704 705 706 707 708
#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
709 710
    IpplTimings::TimerRef efieldTimer_m;        ///< timer for e-field calculation + copy
    IpplTimings::TimerRef averageTimer_m;       ///> timer for average down process
711
    IpplTimings::TimerRef bottomTimer_m;        ///< bottom solver timer
712
    IpplTimings::TimerRef dumpTimer_m;          ///< write SDDS file timer
713 714 715
#endif
};

716 717 718 719 720

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

721
#endif