pfc.cc 3.16 KB
Newer Older
1
2
3
4
5
6
7
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>
#include <ctime>
#include <cmath>

8
9
10
11
12
#include <amdis/AMDiS.hpp>
#include <amdis/AdaptInstationary.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemInstat.hpp>
#include <amdis/Terms.hpp>
13

14
#include "MTLPfcPrecon.hpp"
15
16
17

using namespace AMDiS;

18
19
20
// using Mesh = Dune::YaspGrid<AMDIS_DIM>;
using Mesh = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;

21
using PfcParam = LagrangeBasis<Mesh::LeafGridView, 2, 2, 2>;
22
23


24
25
// 3 component with polynomial degree 1
using PfcProblem = ProblemStat<PfcParam>;
26
using PfcProblemInstat = ProblemInstat<PfcParam>;
27

28
29
30
31
32
using MatrixType = PfcProblem::SystemMatrixType::MultiMatrix;
using VectorType = PfcProblem::SystemVectorType::BaseVector;

using Precon = MTLPfcPrecon<MatrixType,VectorType>;

33
int main(int argc, char** argv)
34
35
{
    AMDiS::init(argc, argv);
36

37
38
    auto creator = new Precon::Creator;
    CreatorMap<typename Precon::PreconBase>::addCreator("pfc", creator);
39

40
41
    PfcProblem prob("pfc");
    prob.initialize(INIT_ALL);
42

43
    PfcProblemInstat probInstat("pfc", prob);
44
    probInstat.initialize(INIT_UH_OLD);
45

46
    AdaptInfo adaptInfo("adapt", prob.getNumComponents());
47

48
49
50
    double r = -0.4;
    double psi_mean = -0.3;
    double M0 = 1.0;
51

52
53
54
    auto& Mu  = prob.getSolution<0>();
    auto& Psi = prob.getSolution<1>();
    auto& Nu  = prob.getSolution<2>();
55

56
    using Op = PfcProblem::ElementOperator;
57

58
    // < [-(1+r) - 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
59
    auto opLhs01_ = Op::zot( -(1.0 + r) );
Praetorius, Simon's avatar
Praetorius, Simon committed
60
61
    opLhs01_.addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
    opLhs01_.addSOT( 2.0 );
62
    auto opLhs01 = std::make_pair(opLhs01_, true);
63

64
    // < -2*psi^3, v >
65
    auto opRhs0 = Op::zot(valueOfFunc(Psi, [](auto psi) { return -2.0 * Math::pow<3>(psi); }, 3));
66

67
    // < psi, v >
68
    auto opRhs1 = Op::zot(valueOf(Psi));
69

70
    // < u, v >
71
    auto opM = std::make_pair(Op::zot(1.0), false);
72

73
    // < grad(u), grad(v) >
74
    auto opL = std::make_pair(Op::sot(1.0), false);
75

76
    // < M0*tau * grad(u), grad(v) >
77
    double* tauPtr = adaptInfo.getTimestepPtr();
78
    auto opL0 = std::make_pair(Op::sot([M0, tauPtr](auto /*x*/) { return M0 * (*tauPtr); }), false);
79
80


81
    // === Define the linear system ===
82

83
84
85
    prob.addMatrixOperator({{ {{0,0}, opM},  {{0,1}, opLhs01}, {{0,2}, opL},
                              {{1,0}, opL0}, {{1,1}, opM},
                                             {{2,1}, opL},     {{2,2}, opM}}});
86
87

    prob.addVectorOperator({{ {0, opRhs0},
88
                              {1, opRhs1} }});
89

90
    // === set the initial solution ===
91
92
    Mu = 0.0;
    Nu = 0.0;
93
    std::srand( 12345 );
94
    Psi.interpol([psi_mean](auto const& x) { return 0.5*(std::rand()/double(RAND_MAX) - 0.5) + psi_mean; });
95

96
    // === configure preconditioner, if necessary ===
97
    if (creator->precon) {
98
      creator->precon->setData(tauPtr, M0);
99
    }
100

101
102
    AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
    adapt.adapt();
103

104
105
106
    msg("|mu|  = ", two_norm(Mu.getVector()), ", ",
        "|nu|  = ", two_norm(Nu.getVector()), ", ",
        "|psi| = ", two_norm(Psi.getVector()));
107

108
109
110
    AMDiS::finalize();
    return 0;
}