AmrSmoother.cpp 3.27 KB
Newer Older
1 2 3 4 5 6 7 8
// Source file of the AmrSmoother class,
//   interface to Ifpack2 smoothers.
//
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// OPAL is licensed under GNU GPL version 3.

9
#include "AmrSmoother.h"
10 11

#include <map>
12 13
#include <utility>

14 15 16
#include "Utilities/OpalException.h"
#include "Utilities/Util.h"

17 18
AmrSmoother::AmrSmoother(const Teuchos::RCP<const matrix_t>& A,
                         const Smoother& smoother,
19
                         lo_t nSweeps)
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
{
    const std::string type = "RELAXATION";
    
    Ifpack2::Factory factory;
    prec_mp = factory.create(type, A);
    
    params_mp = Teuchos::rcp( new Teuchos::ParameterList );
    
    this->initParameter_m(smoother, nSweeps);
    
    
    prec_mp->setParameters(*params_mp);
    prec_mp->initialize();
    prec_mp->compute();
}


AmrSmoother::~AmrSmoother() {
    prec_mp = Teuchos::null;
    params_mp = Teuchos::null;
}


void AmrSmoother::smooth(const Teuchos::RCP<vector_t>& x,
                         const Teuchos::RCP<matrix_t>& A,
                         const Teuchos::RCP<vector_t>& b)
{
    prec_mp->apply(*b, *x, Teuchos::NO_TRANS,
                   Teuchos::ScalarTraits<scalar_t>::one(),
                   Teuchos::ScalarTraits<scalar_t>::zero());
}


53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
AmrSmoother::Smoother
AmrSmoother::convertToEnumSmoother(const std::string& smoother) {
    std::map<std::string, Smoother> map;
    
    map["GS"]     = Smoother::GAUSS_SEIDEL;
    map["SGS"]    = Smoother::SGS;
    map["JACOBI"] = Smoother::JACOBI;
    
    auto sm = map.find(Util::toUpper(smoother));
    
    if ( sm == map.end() )
        throw OpalException("AmrMultiGrid::convertToEnumNorm_m()",
                            "No smoother '" + smoother + "'.");
    return sm->second;
}


70
void AmrSmoother::initParameter_m(const Smoother& smoother,
71
                                  lo_t nSweeps)
72 73 74 75 76 77
{
    if ( params_mp == Teuchos::null )
        params_mp = Teuchos::rcp( new Teuchos::ParameterList );
    
    
    std::string type = "";
78 79
    scalar_t damping = 1.0;
    std::pair<bool, scalar_t> l1 = std::make_pair(true, 1.5);
80 81
    
    bool backward = false;
82
    std::pair<bool, scalar_t> fix = std::make_pair(true, 1.0e-5);
83 84 85 86 87 88 89 90 91 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
    bool check = true;
    
    switch ( smoother ) {
        
        case GAUSS_SEIDEL:
        {
            type = "Gauss-Seidel";
            backward = false;
            damping = 1.0;
            break;
        }
        case SGS:
        {
            type = "Symmetric Gauss-Seidel";
            damping = 1.0;
            break;
        }
        case JACOBI:
        {
            type = "Jacobi";
            damping = 6.0 / 7.0;
            break;
        }
        default:
            break;
    };
    
    params_mp->set("relaxation: type", type);
    params_mp->set("relaxation: sweeps", nSweeps);
    params_mp->set("relaxation: zero starting solution", false);
    params_mp->set("relaxation: damping factor", damping);
    params_mp->set("relaxation: use l1", l1.first);
    params_mp->set("relaxation: l1 eta", l1.second);
    params_mp->set("relaxation: backward mode", backward);
    params_mp->set("relaxation: fix tiny diagonal entries", fix.first);
    params_mp->set("relaxation: min diagonal value", fix.second);
    params_mp->set("relaxation: check diagonal entries", check);
}