heat.cc 5.92 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
class Heat : public ProblemInstat
47
48
{
public:
49
50
  Heat(ProblemStat &heatSpace)
    : ProblemInstat("heat", heatSpace)
51
52
53
  {
    // init theta scheme
    theta = -1.0;
54
    Parameters::get(name + "->theta", theta);
55
56
57
58
59
60
    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
    ProblemInstat::setTime(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
    rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
  }
72

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

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

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

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

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

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

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

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

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

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

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

private:
127
  /// Used for theta scheme.
128
129
  double theta;

130
  /// theta - 1
131
132
  double theta1;

133
  /// time for right hand side functions.
134
135
  double rhsTime;

136
  /// Pointer to boundary function. Needed for initial problem.
137
138
139
140
141
142
143
144
145
  AbstractFunction<double, WorldVector<double> > *exactSolution;
};

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

int main(int argc, char** argv)
{
146
  FUNCNAME("main");
147

148
  AMDiS::init(argc, argv);
149

150
  // ===== create and init stationary problem =====
151
  ProblemStat heatSpace("heat->space");
152
153
  heatSpace.initialize(INIT_ALL);

154
155

  // ===== create instationary problem =====
156
157
  Heat heat(heatSpace);
  heat.initialize(INIT_ALL);
158
159

  // create adapt info
160
  AdaptInfo adaptInfo("heat->adapt", heatSpace.getNumComponents());
161
162

  // create initial adapt info
163
  AdaptInfo adaptInfoInitial("heat->initial->adapt");
164
165

  // create instationary adapt
166
167
168
169
170
171
  AdaptInstationary adaptInstat("heat->adapt",
				heatSpace,
				adaptInfo,
				heat,
				adaptInfoInitial);

172
173

  // ===== create rhs functions =====
174
  int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
175
  F *rhsFct = new F(degree);
176
  rhsFct->setTimePtr(heat.getRhsTimePtr());
177

178
179
180
181
182
183

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

  // create laplace
184
  Operator A(heatSpace.getFeSpace());
185
186
  A.addSecondOrderTerm(new Simple_SOT);
  A.setUhOld(heat.getOldSolution(0));
187
  if (*(heat.getThetaPtr()) != 0.0)
188
    heatSpace.addMatrixOperator(A, 0, 0, heat.getThetaPtr(), &one);
189
  if (*(heat.getTheta1Ptr()) != 0.0)
190
    heatSpace.addVectorOperator(A, 0, heat.getTheta1Ptr(), &zero);
191
192

  // create zero order operator
193
  Operator C(heatSpace.getFeSpace());
194
  C.addZeroOrderTerm(new Simple_ZOT);
195
196
197
  C.setUhOld(heat.getOldSolution(0));
  heatSpace.addMatrixOperator(C, 0, 0, heat.getInvTau(), heat.getInvTau());
  heatSpace.addVectorOperator(C, 0, heat.getInvTau(), heat.getInvTau());
198
199

  // create RHS operator
200
  Operator F(heatSpace.getFeSpace());
201
  F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
202
  heatSpace.addVectorOperator(F, 0);
203
204


205
206
207
208
  // ===== create boundary functions =====
  G *boundaryFct = new G;
  boundaryFct->setTimePtr(heat.getTime());
  heat.setExactSolution(boundaryFct);
209
  DOFVector<double> G_dof(heatSpace.getFeSpace(), "G");
Praetorius, Simon's avatar
Praetorius, Simon committed
210
211
  G_dof.interpol(boundaryFct);
  heatSpace.addDirichletBC(1, 0, 0, &G_dof);
212
213


214
  // ===== start adaption loop =====
215
  int errorCode = adaptInstat.adapt();
216

217
218
  AMDiS::finalize();

219
220
  return errorCode;
}