heat.cc 7.04 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

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

/** \brief
 * Dirichlet boundary function
 */
class G : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
  MEMORY_MANAGED(G);

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
22
23
  double operator()(const WorldVector<double>& x) const {
    return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x));
Thomas Witkowski's avatar
Thomas Witkowski committed
24
  }
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
};

/** \brief
 * RHS function
 */
class F : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
  MEMORY_MANAGED(F);

  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
41
  double operator()(const WorldVector<double>& x) const {
42
    int dim = x.getSize();
43
44
45
46
    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
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
};

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

/** \brief
 * Instationary problem
 */
class Heat : public ProblemInstatScal
{
public:
  MEMORY_MANAGED(Heat);

  /** \brief
   *  Constructor
   */
  Heat(ProblemScal *heatSpace)
    : 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);
    theta1 = theta - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
78
  }
79
80
81
82
83
84
85
86


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

  /** \brief
   * set the time in all needed functions!
   */
  void setTime(AdaptInfo *adaptInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
87
    rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
88
89
    boundaryTime = adaptInfo->getTime();
    tau1 = 1.0 / adaptInfo->getTimestep();    
Thomas Witkowski's avatar
Thomas Witkowski committed
90
  }
91
92
93
94

  void closeTimestep(AdaptInfo *adaptInfo) {
    ProblemInstatScal::closeTimestep(adaptInfo);
    WAIT;
Thomas Witkowski's avatar
Thomas Witkowski committed
95
  }
96
97
98
99
100
101

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

  /** \brief
   * Used by \ref problemInitial to solve the system of the initial problem
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
102
  void solve(AdaptInfo *adaptInfo, bool) {
103
    problemStat->getSolution()->interpol(exactSolution);
Thomas Witkowski's avatar
Thomas Witkowski committed
104
  }
105
106
107
108
109

  /** \brief
   * Used by \ref problemInitial to do error estimation for the initial
   * problem.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
110
  void estimate(AdaptInfo *adaptInfo) {
111
112
113
114
115
116
117
    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
118
  }
119
120
121
122
123
124

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

  /** \brief
   * Sets \ref exactSolution;
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
125
  void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct) {
126
127
128
129
130
131
132
133
    exactSolution = fct;
  } 

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

  /** \brief
   * Returns pointer to \ref theta.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
134
  double *getThetaPtr() { 
135
    return &theta; 
Thomas Witkowski's avatar
Thomas Witkowski committed
136
  }
137
138
139
140

  /** \brief
   * Returns pointer to \ref theta1.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
141
  double *getTheta1Ptr() { 
142
    return &theta1; 
Thomas Witkowski's avatar
Thomas Witkowski committed
143
  }
144
145
146
147

  /** \brief
   * Returns pointer to \ref tau1
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
148
  double *getTau1Ptr() { 
149
    return &tau1; 
Thomas Witkowski's avatar
Thomas Witkowski committed
150
  }
151
152
153
154

  /** \brief
   * Returns pointer to \ref rhsTime.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
155
  double *getRHSTimePtr() { 
156
    return &rhsTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
157
  }
158
159
160
161

  /** \brief
   * Returns pointer to \ref theta1.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
162
  double *getBoundaryTimePtr() {
163
    return &boundaryTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
164
  }
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288

private:
  /** \brief
   * Used for theta scheme.
   */
  double theta;

  /** \brief
   * theta - 1
   */
  double theta1;

  /** \brief
   * 1.0 / timestep
   */
  double tau1;

  /** \brief
   * time for right hand side functions.
   */
  double rhsTime;

  /** \brief
   * time for boundary functions.
   */
  double boundaryTime;

  /** \brief
   * Pointer to boundary function. Needed for initial problem.
   */
  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);

  // ===== create and init stationary problem =====
  ProblemScal *heatSpace = NEW ProblemScal("heat->space");
  heatSpace->initialize(INIT_ALL);

  // ===== create instationary problem =====
  Heat *heat = new Heat(heatSpace);
  heat->initialize(INIT_ALL);

  // create adapt info
  AdaptInfo *adaptInfo = NEW AdaptInfo("heat->adapt");

  // create initial adapt info
  AdaptInfo *adaptInfoInitial = NEW AdaptInfo("heat->initial->adapt");

  // create instationary adapt
  AdaptInstationary *adaptInstat =
    NEW AdaptInstationary("heat->adapt",
			  heatSpace,
			  adaptInfo,
			  heat,
			  adaptInfoInitial);

  // ===== create boundary functions =====
  G *boundaryFct = NEW G;
  boundaryFct->setTimePtr(heat->getBoundaryTimePtr());
  heat->setExactSolution(boundaryFct);

  heatSpace->addDirichletBC(DIRICHLET, boundaryFct);

  // ===== create rhs functions =====
  int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree();
  F *rhsFct = NEW F(degree);
  rhsFct->setTimePtr(heat->getRHSTimePtr());

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

  // create laplace
  Operator *A = NEW Operator(Operator::MATRIX_OPERATOR | 
			     Operator::VECTOR_OPERATOR,
			     heatSpace->getFESpace());

  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);

  // create zero order operator
  Operator *C = NEW Operator(Operator::MATRIX_OPERATOR | 
			     Operator::VECTOR_OPERATOR,
			     heatSpace->getFESpace());

  C->addZeroOrderTerm(NEW Simple_ZOT);

  C->setUhOld(heat->getOldSolution());

  heatSpace->addMatrixOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
  heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());

  // create RHS operator
  Operator *F = NEW Operator(Operator::VECTOR_OPERATOR,
			     heatSpace->getFESpace());

  F->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct));

  heatSpace->addVectorOperator(F);

  // ===== start adaption loop =====
  int errorCode = adaptInstat->adapt();

  return errorCode;
}