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

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <ctime>
#include <cmath>

11
#include <dune/amdis/AdaptInstationary.hpp>
12
13
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
14
#include <dune/amdis/ProblemInstat.hpp>
15
16
17
18
19
20
21
22

#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif

23
// #include "PfcPrecon.hpp"
24
25
26

using namespace AMDiS;

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

class PfcParam
{
  template <class Mesh, int deg>
  using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
  
  template <class Mesh, int... degs>
  using FeSpaceTuple = std::tuple<Lagrange<Mesh, degs>...>;
  
public:
  static constexpr int dim = DIM;
  static constexpr int dimworld = DOW;
  static constexpr int nComponents = 3;
  
  // default values
  using Mesh = Dune::YaspGrid<dim>;
  using FeSpaces = FeSpaceTuple<Mesh, 2, 2, 2>;
};

46
// 3 component with polynomial degree 1
47
// using PfcParam   = ProblemStatTraits<DIM, DOW, 1, 1, 1>;
48
using PfcProblem = ProblemStat<PfcParam>;
49
using PfcProblemInstat = ProblemInstat<PfcParam>;
50

51
int main(int argc, char** argv)
52
53
54
55
56
57
{
    AMDiS::init(argc, argv);
    
    PfcProblem prob("pfc");
    prob.initialize(INIT_ALL);
    
58
    PfcProblemInstat probInstat("pfc", prob);
59
    probInstat.initialize(INIT_UH_OLD);
60
    
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    AdaptInfo adaptInfo("adapt");
    
    double r = -0.4;
    double psi_mean = -0.3;
    double M0 = 1.0;
    
    auto& Mu  = prob.getSolution<0>();
    auto& Psi = prob.getSolution<1>();
    auto& Nu  = prob.getSolution<2>();
    
    using Op = PfcProblem::OperatorType;
    Op opLhs00, opLhs01, opLhs02; 
    opLhs00.addZOT( constant(1.0) );
    opLhs01.addZOT( constant(-(1.0 + r)) );
    opLhs01.addZOT( valueOfFunc(Psi, [](auto psi) { return 3.0 * std::pow(psi,2); }, 2) );
    opLhs01.addSOT( constant(2.0) );
    opLhs02.addSOT( constant(1.0) );
    
    Op opLhs10, opLhs11; 
80
    opLhs10.addSOT( constant(M0) );
81
82
83
84
85
86
87
88
89
90
91
92
93
    opLhs11.addZOT( constant(1.0) );
    
    Op opLhs21, opLhs22; 
    opLhs21.addSOT( constant(1.0) );
    opLhs22.addZOT( constant(1.0) );
    
    Op opRhs0, opRhs1;  
    opRhs0.addZOT( valueOfFunc(Psi, [](auto psi) { return -2.0 * std::pow(psi,3); }, 3) );
    opRhs1.addZOT( valueOf(Psi) );
    
    prob.addMatrixOperator(opLhs00, 0, 0);
    prob.addMatrixOperator(opLhs01, 0, 1);
    prob.addMatrixOperator(opLhs02, 0, 2);
94
    prob.addMatrixOperator(opLhs10, 1, 0, adaptInfo.getTimestepPtr());
95
96
97
    prob.addMatrixOperator(opLhs11, 1, 1);
    prob.addMatrixOperator(opLhs21, 2, 1);
    prob.addMatrixOperator(opLhs22, 2, 2);
98
    
99
100
101
102
103
104
105
    prob.addVectorOperator(opRhs0, 0);
    prob.addVectorOperator(opRhs1, 1);
    
    // initial solution
    Mu.getVector() = 0.0;
    Nu.getVector() = 0.0;
    std::srand( time(0) );
106
    Psi.interpol([psi_mean](auto const& x) { return (std::rand()/double(RAND_MAX) - 0.5) + psi_mean; });
107
    
108
109
//     using PreconType = PfcPrecon<typename PfcProblem::SystemMatrixType::MultiMatrix, typename PfcProblem::SystemVectorType::MultiVector>;
//     PreconType precon(prob.getSystemMatrix().getMatrix(), &tau, M0);
110
    
111
112
113
    AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
    adapt.adapt();
    
114
115
116
    AMDiS::finalize();
    return 0;
}