RosenbrockStationary.cc 5.22 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
#include "time/RosenbrockStationary.h"
23
#include "io/VtkWriter.h"
24
#include "ProblemStat.h"
25
#include "SystemVector.h"
26
// #include "solver/OEMSolver.h"
27
#include "Debug.h"
28

29
30
31
32
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MeshDistributor.h"
#endif

33
34
35
36
37
38
39
40
namespace AMDiS {

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

41

42
  void RosenbrockStationary::init()
43
44
45
  {
    stageSolution = new SystemVector(*solution);
    unVec = new SystemVector(*solution);
46
    timeRhsVec = new SystemVector(*solution);
47
48
    newUn = new SystemVector(*solution);
    tmp = new SystemVector(*solution);
49
    lowSol = new SystemVector(*solution);    
50
51
52

    stageSolution->set(0.0);
    unVec->set(0.0);
53
    
54
    stageSolutions.resize(rm->getStages());
55
    for (int i = 0; i < rm->getStages(); i++) {
56
      stageSolutions[i] = new SystemVector(*solution);
57
58
      stageSolutions[i]->set(0.0);
    }
59
60
61
62
63
64
  }


  void RosenbrockStationary::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
					       bool asmMatrix, bool asmVector)
  {        
65
66
67
68
    FUNCNAME("RosenbrockStationary::buildAfterCoarsen()");

    TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n");

69
70
    if (first) {
      first = false;
71
      *unVec = *solution;      
72
73
74
75
    }
    
    *newUn = *unVec;    
    *lowSol = *unVec;
76
77
78
79
for (int i = 0; i < nComponents; i++) {
  std::cout <<"RRDDRX "<< i<<" "<< lowSol->getDOFVector(i)->Int() <<" "<< 
    newUn->getDOFVector(i)->Int()<<" "<< *tauPtr<<"\n";
 }
80
    for (int i = 0; i < rm->getStages(); i++) {      
81
82
83
      *stageSolution = *unVec;
      for (int j = 0; j < i; j++) {
	*tmp = *(stageSolutions[j]);
84
	*tmp *= rm->getA(i, j);
85
86
	*stageSolution += *tmp;
      }
87
88
89
90
91
92

      for (unsigned int j = 0; j < boundaries.size(); j++) {
	boundaries[j].vec->interpol(boundaries[j].fct);
	*(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].row));
      }

93
      timeRhsVec->set(0.0);
94
      for (int j = 0; j < i; j++) {
95
96
97
	*tmp = *(stageSolutions[j]);
	*tmp *= (rm->getC(i, j) / *tauPtr);
	*timeRhsVec += *tmp;
98
      }
99

100
      ProblemStat::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
101
102
103
104
105
#if defined HAVE_PARALLEL_PETSC
      // TODO: Problems with reuse of Matrix with parallel PETSC-Solvers
      //       Thus, Rosenbrock not efficient but working (Rainer)
      ProblemStat::solve(adaptInfo, true , false);
#else
106
      ProblemStat::solve(adaptInfo, i == 0, i + 1 < rm->getStages());
107
108
#endif
  
109

110
111
112
      *(stageSolutions[i]) = *solution;
      
      *tmp = *solution;
113
      *tmp *= rm->getM1(i);
114

115
      *newUn += *tmp;
116

117
      *tmp = *solution;
118
      *tmp *= rm->getM2(i);
119
120
121
122
123
      *lowSol += *tmp;
    }
    
    for (int i = 0; i < nComponents; i++) {
      (*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
124
      adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->L2Norm(), i);
125
    }   
126
127
128
  }


129
  void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool, bool)
130
131
  {}

132
133
134
135
136
137

  void RosenbrockStationary::addOperator(Operator &op, int row, int col, 
					 double *factor, double *estFactor)
  {
    FUNCNAME("RosenbrockStationary::addOperator()");

138
    TEST_EXIT(op.getUhOld() == nullptr)("UhOld is not allowed to be set!\n");
139
140

    op.setUhOld(stageSolution->getDOFVector(col));
Thomas Witkowski's avatar
Thomas Witkowski committed
141
    ProblemStat::addVectorOperator(op, row, factor, estFactor);
142
143
144
145
146
147
148
149
  }
  

  void RosenbrockStationary::addJacobianOperator(Operator &op, int row, int col, 
						 double *factor, double *estFactor)
  {
    FUNCNAME("RosenbrockStationary::addJacobianOperator()");
    
150
151
    TEST_EXIT(factor == nullptr)("Not yet implemented!\n");
    TEST_EXIT(estFactor == nullptr)("Not yet implemented!\n");
152

153
    ProblemStat::addMatrixOperator(op, row, col, &minusOne, &minusOne);
154
155
156
157
158
159
160
  }


  void RosenbrockStationary::addTimeOperator(int row, int col)
  {
    FUNCNAME("RosenbrockStationary::addTimeOperator()");

Thomas Witkowski's avatar
Thomas Witkowski committed
161
    TEST_EXIT(invTauGamma)("This should not happen!\n");
162
163
164

    Operator *op = new Operator(componentSpaces[row], componentSpaces[col]);
    op->addZeroOrderTerm(new Simple_ZOT);
Thomas Witkowski's avatar
Thomas Witkowski committed
165
    ProblemStat::addMatrixOperator(op, row, col, invTauGamma, invTauGamma);
166
167

    Operator *opRhs = new Operator(componentSpaces[row]);
Thomas Witkowski's avatar
Thomas Witkowski committed
168
    opRhs->addZeroOrderTerm(new VecAtQP_ZOT(timeRhsVec->getDOFVector(col)));
169
    ProblemStat::addVectorOperator(opRhs, row);
170
171
  }

172
173
174
175
176
177
178
179
  
  void RosenbrockStationary::addDirichletBC(BoundaryType type, int row, int col,
					    AbstractFunction<double, WorldVector<double> > *fct)
  {
    DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
    RosenbrockBoundary bound = {fct, vec, row, col};
    boundaries.push_back(bound);

180
    ProblemStat::addDirichletBC(type, row, col, vec);
181
182
  }

183
}