RosenbrockStationary.cc 2.91 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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
121
#include "time/RosenbrockStationary.h"
#include "ProblemVec.h"
#include "SystemVector.h"
#include "OEMSolver.h"

namespace AMDiS {

  void RosenbrockStationary::acceptTimestep()
  {
    *solution = *newUn;
    *unVec = *newUn;
  }

  void RosenbrockStationary::initRB()
  {
    stages = 3;
    
    stageSolution = new SystemVector(*solution);
    unVec = new SystemVector(*solution);
    newUn = new SystemVector(*solution);
    tmp = new SystemVector(*solution);
    lowSol = new SystemVector(*solution);
    
    stageSolutions.resize(stages);
    for (int i = 0; i < stages; i++)
      stageSolutions[i] = new SystemVector(*solution);
    
    phiSum = new DOFVector<double>(feSpaces[0], "phiSum");
    tmpDof = new DOFVector<double>(feSpaces[0], "phiSum");
    
    a.resize(stages);
    for (int i = 0; i < stages; i++) {
      a[i].resize(stages);
      for (int j = 0; j < stages; j++)
	a[i][j] = 0.0;
    }
    a[0][0] = 0.0;
    a[1][0] = 1.605996252195329e+00;
    a[1][1] = 0.0;
    a[2][0] = 1.605996252195329e+00;
    a[2][1] = 0.0;
    a[2][2] = 0.0;
    
    c.resize(stages);
    for (int i = 0; i < stages; i++) {
      c[i].resize(stages);
      for (int j = 0; j < stages; j++)
	c[i][j] = 0.0;
    }
    
    c[0][0] =  2.294280360279042e+00;
    c[1][0] = -8.874044410657833e-01;
    c[1][1] =  2.294280360279042e+00;
    c[2][0] = -2.398747971635036e+01;
    c[2][1] = -5.263722371562129e+00;
    c[2][2] =  2.294280360279042e+00;
    
    m.resize(stages);
    m[0] =  2.236727045296590e+00;
    m[1] =  2.250067730969644e+00;
    m[2] = -2.092514044390320e-01;
    
    m2.resize(stages);
    m2[0] = 2.059356167645940e+00;
    m2[1] = 1.694014319346528e-01;
    m2[2] = 0.0;
  }


  void RosenbrockStationary::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
					       bool asmMatrix, bool asmVector)
  {        
    if (first) {
      first = false;
      *unVec = *solution;
    }
    
    *newUn = *unVec;    
    *lowSol = *unVec;
    
    for (int i = 0; i < stages; i++) {      
      *stageSolution = *unVec;
      for (int j = 0; j < i; j++) {
	*tmp = *(stageSolutions[j]);
	*tmp *= a[i][j];
	*stageSolution += *tmp;
      }
      
      phiSum->set(0.0);
      for (int j = 0; j < i; j++) {
	*tmpDof = *(stageSolutions[j]->getDOFVector(0));
	*tmpDof *= (c[i][j] / *tauPtr);
	*phiSum += *tmpDof;
      }
      
      ProblemVec::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
      solver->setMultipleRhs(i != 0);
      ProblemVec::solve(adaptInfo);
      
      *(stageSolutions[i]) = *solution;
      
      *tmp = *solution;
      *tmp *= m[i];
      *newUn += *tmp;
      
      *tmp = *solution;
      *tmp *= m2[i];
      *lowSol += *tmp;
    }
    
    for (int i = 0; i < nComponents; i++) {
      (*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
      adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i);
    }
  }


  void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
  {}

}