heat.cc 6.14 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

10
/// Dirichlet boundary function
11
12
13
14
class G : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
15
16
17
18
  
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
19
    return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x));
Thomas Witkowski's avatar
Thomas Witkowski committed
20
  }
21
22
};

23
/// RHS function
24
25
26
27
class F : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
28
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
29

30
31
32
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
33
    int dim = Global::getGeo(WORLD);
34
35
36
37
    double r2 = x * x;
    double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2);
    double ut = M_PI * cos(M_PI * (*timePtr)) * exp(-10.0 * r2);
    return ut -(400.0 * r2 - 20.0 * dim) * ux;
Thomas Witkowski's avatar
Thomas Witkowski committed
38
  }
39
40
41
42
43
44
};

// ===========================================================================
// ===== instationary problem ================================================
// ===========================================================================

45
/// Instationary problem
46
47
48
class Heat : public ProblemInstatScal
{
public:
49
  Heat(ProblemScal &heatSpace)
50
51
52
53
54
55
56
57
58
59
60
    : ProblemInstatScal("heat", heatSpace)
  {
    // init theta scheme
    theta = -1.0;
    GET_PARAMETER(0, name + "->theta", "%f", &theta);
    TEST_EXIT(theta >= 0)("theta not set!\n");
    if (theta == 0.0) {
      WARNING("You are using the explicit Euler scheme\n");
      WARNING("Use a sufficiently small time step size!!!\n");
    }
    MSG("theta = %f\n", theta);
61
    theta1 = theta - 1.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
62
  }
63
64
65

  // ===== ProblemInstatBase methods ===================================

66
67
68
  /// set the time in all needed functions!
  void setTime(AdaptInfo *adaptInfo) 
  {
69
70
    ProblemInstat::setTime(adaptInfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
71
    rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
72
    boundaryTime = adaptInfo->getTime();
Thomas Witkowski's avatar
Thomas Witkowski committed
73
  }
74

75
76
  void closeTimestep(AdaptInfo *adaptInfo) 
  {
77
78
    ProblemInstatScal::closeTimestep(adaptInfo);
    WAIT;
Thomas Witkowski's avatar
Thomas Witkowski committed
79
  }
80
81
82

  // ===== initial problem methods =====================================

83
84
85
  /// Used by \ref problemInitial to solve the system of the initial problem
  void solve(AdaptInfo *adaptInfo) 
  {
86
    problemStat->getSolution()->interpol(exactSolution);
Thomas Witkowski's avatar
Thomas Witkowski committed
87
  }
88

89
90
91
  /// Used by \ref problemInitial to do error estimation for the initial problem.
  void estimate(AdaptInfo *adaptInfo) 
  {
92
93
94
95
96
97
    double errMax, errSum;

    errSum = Error<double>::L2Err(*exactSolution,
				  *(problemStat->getSolution()), 0, &errMax, false);
    adaptInfo->setEstSum(errSum, 0);
    adaptInfo->setEstMax(errMax, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
98
  }
99
100
101

  // ===== setting methods ===============================================

102
103
104
  /// Sets \ref exactSolution;
  void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
105
106
107
108
109
    exactSolution = fct;
  } 

  // ===== getting methods ===============================================

110
111
112
  /// Returns pointer to \ref theta.
  double *getThetaPtr() 
  { 
113
    return &theta; 
Thomas Witkowski's avatar
Thomas Witkowski committed
114
  }
115

116
117
118
  /// Returns pointer to \ref theta1.
  double *getTheta1Ptr() 
  { 
119
    return &theta1; 
Thomas Witkowski's avatar
Thomas Witkowski committed
120
  }
121

122
123
124
  /// Returns pointer to \ref rhsTime.
  double *getRHSTimePtr() 
  { 
125
    return &rhsTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
126
  }
127

128
129
130
  /// Returns pointer to \ref theta1.
  double *getBoundaryTimePtr() 
  {
131
    return &boundaryTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
132
  }
133
134

private:
135
  /// Used for theta scheme.
136
137
  double theta;

138
  /// theta - 1
139
140
  double theta1;

141
  /// time for right hand side functions.
142
143
  double rhsTime;

144
  /// time for boundary functions.
145
146
  double boundaryTime;

147
  /// Pointer to boundary function. Needed for initial problem.
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
  AbstractFunction<double, WorldVector<double> > *exactSolution;
};

// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

int main(int argc, char** argv)
{
  // ===== check for init file =====
  TEST_EXIT(argc > 1)("usage: heat initfile\n");

  // ===== init parameters =====
  Parameters::init(false, argv[1]);
  Parameters::readArgv(argc, argv);

164

165
  // ===== create and init stationary problem =====
166
167
168
  ProblemScal heatSpace("heat->space");
  heatSpace.initialize(INIT_ALL);

169
170

  // ===== create instationary problem =====
171
172
  Heat heat(heatSpace);
  heat.initialize(INIT_ALL);
173
174

  // create adapt info
175
  AdaptInfo adaptInfo("heat->adapt");
176
177

  // create initial adapt info
178
  AdaptInfo adaptInfoInitial("heat->initial->adapt");
179
180

  // create instationary adapt
181
182
183
184
185
186
  AdaptInstationary adaptInstat("heat->adapt",
				heatSpace,
				adaptInfo,
				heat,
				adaptInfoInitial);

187

188
189
190
191
192
193
194
  // ===== create boundary functions =====
  G *boundaryFct = new G;
  boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
  heat.setExactSolution(boundaryFct);
  heatSpace.addDirichletBC(1, boundaryFct);


195
  // ===== create rhs functions =====
196
  int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
197
  F *rhsFct = new F(degree);
198
199
  rhsFct->setTimePtr(heat.getRHSTimePtr());

200
201
202
203
204
205

  // ===== create operators =====
  double one = 1.0;
  double zero = 0.0;

  // create laplace
206
  Operator A(heatSpace.getFeSpace());
207
208
209
210
211
212
  A.addSecondOrderTerm(new Laplace_SOT);
  A.setUhOld(heat.getOldSolution());
  if (*(heat.getThetaPtr()) != 0.0)
    heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one);
  if (*(heat.getTheta1Ptr()) != 0.0)
    heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);
213
214

  // create zero order operator
215
  Operator C(heatSpace.getFeSpace());
216
217
  C.addZeroOrderTerm(new Simple_ZOT);
  C.setUhOld(heat.getOldSolution());
218
219
  heatSpace.addMatrixOperator(C, heat.getInvTau(), heat.getInvTau());
  heatSpace.addVectorOperator(C, heat.getInvTau(), heat.getInvTau());
220
221

  // create RHS operator
222
  Operator F(heatSpace.getFeSpace());
223
224
  F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
  heatSpace.addVectorOperator(F);
225
226
227


  // ===== start adaption loop =====
228
  int errorCode = adaptInstat.adapt();
229
230
231

  return errorCode;
}