couple.tex 12 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
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
\section{Coupled problems}
\label{s:couple}
In this example, we solve the same problem as in Section \ref{s:vecellipt}, but here we treat the two equations as two coupled problems. The main difference is that the equations now aren't assembled into the same large system of equations, but into two separated systems of equations, that have to be solved separately. We define the two problems
\begin{eqnarray}
-\Delta u &=& f ~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
u &=& g ~~~ \mbox{on } \partial\Omega
\end{eqnarray}
and
\begin{equation}
v = u.
\end{equation}
We first solve the first problem and then use its solution to solve the second problem. This happens in every iteration of the adaptation loop. Both problems should use the same mesh. Mesh adaptation is done by the first problem. So one iteration now looks like illustrated in Figure \ref{f:coupled iteration}.  
\begin{figure}
\center
\includegraphics[width=0.9\textwidth]{fig/coupled_iteration.pdf}
\caption{State diagram of the coupled iteration.}
\label{f:coupled iteration}
\end{figure} 

\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
\hline 
& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
\hline 
{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt couple.cc} \\
\hline 
{\bf parameter file} & \tt init/ & \tt couple.dat.1d & \tt couple.dat.2d & \tt couple.dat.3d \\
\hline 
{\bf macro file} & \tt macro/ & \tt macro.stand.1d & \tt macro.stand.2d & \tt macro.stand.3d  \\
\hline 
{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt couple.mesh, couple.dat} \\
\hline 
\end{tabular}
\caption{Files of the {\tt couple} example.}
\end{table}

\subsection{Source code}
In the previous examples, the iteration was implemented always by the corresponding problem class. In this example, we want to couple two problems within one iteration, so the default implementation can't be used. For that reason, we define our own coupled iteration class\\ \verb+MyCoupledIteration+ which implements the \verb+ProblemIterationInterface+:

\begin{lstlisting}{}
class MyCoupledIteration : public ProblemIterationInterface
{
public:
  MyCoupledIteration(ProblemStatBase *prob1,
		     ProblemStatBase *prob2)
    : problem1(prob1),
      problem2(prob2)
49
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
52
53
54
55
56
57
58
59
60
61
62
\end{lstlisting}

In the constructor pointers to the two problems are assigned to the private members \verb+problem1+ and \verb+problem2+. Note that the pointers point to the interface \verb+ProblemStatBase+ and not to\\ \verb+ProblemScal+. This leads to a more general implementation. If e.g. two vector valued problems should be coupled in the future, we could use our iteration class without modifications. 

Now, we implement the needed interface methods:  

\begin{lstlisting}{}
  void beginIteration(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("StandardProblemIteration::beginIteration()");
    MSG("\n");
    MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
63
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
64
65
66
67
68
69

  void endIteration(AdaptInfo *adaptInfo) {
    FUNCNAME("StandardProblemIteration::endIteration()");
    MSG("\n");
    MSG("end of iteration %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
70
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
71
72
73
74
75
76
77
78
79
80
\end{lstlisting}

These two functions are called at the beginning and at the end of each iteration. Here, we only prompt some output. 

The method \verb+oneIteration+ is the crucial part of our implementation:

\begin{lstlisting}{}
  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) 
  {
    Flag flag, markFlag;
81
82
    if (toDo.isSet(MARK)) markFlag = problem1->markElements(adaptInfo);
    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
Thomas Witkowski's avatar
Thomas Witkowski committed
83
      flag = problem1->refineMesh(adaptInfo);
84
85
86
    
    if (toDo.isSet(BUILD)) problem1->buildAfterCoarsen(adaptInfo, markFlag);
    if (toDo.isSet(SOLVE)) problem1->solve(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
87
 
88
89
    if (toDo.isSet(BUILD)) problem2->buildAfterCoarsen(adaptInfo, markFlag);
    if (toDo.isSet(SOLVE)) problem2->solve(adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
90

91
    if (toDo.isSet(ESTIMATE)) problem1->estimate(adaptInfo);    
Thomas Witkowski's avatar
Thomas Witkowski committed
92
    return flag;
93
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
\end{lstlisting}

The \verb+toDo+ flag is used by the adaptation loop to determine which parts of the iteration should be performed. The first iteration is always an iteration without mesh adaptation (see Figure \ref{f:stationary loop}). So we start our iteration by marking and adapting the mesh. The mesh and its adaptation is managed by the first problem. So we call \verb+markElements+ and \verb+refineMesh+ of \verb+problem1+. Note that no mesh coarsenings have to be performed in our example. Afterwards, \verb+problem1+ assembles its system of equations by \verb+buildAfterCoarsen+. Assemblage and mesh adaptation are nested operations in AMDiS (\verb+buildBeforeRefine+, \verb+refineMesh+, \verb+buildBeforeCoarsen+,\\\verb+coarsenMesh+,\verb+buildAfterCoarsen+). Here, we implement a simplified version.

After \verb+problem1+ has solved its system of equations, \verb+problem2+ can assemble and solve its equations system using the solution of the first problem as right hand side. In the method \verb+oneIteration+, only the order of method calls is determined. The dependency to the solution of the first problem is created later when the operator for the right hand side of \verb+problem2+ is created. 

After also the second problem computed its solution, \verb+problem1+ does the error estimation (remember: mesh adaptation is managed by \verb+problem1+).
\begin{figure}
\center
\includegraphics[width=0.9\textwidth]{fig/stationary_loop.pdf}
\caption{Stationary adaptation loop.}
\label{f:stationary loop}
\end{figure} 

Now, the access to the coupled problems is implemented and the member variables are defined:

\begin{lstlisting}{}
  int getNumProblems() 
  {
    return 2;
114
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
117
118

  ProblemStatBase *getProblem(int number = 0) 
  {
    FUNCNAME("CoupledIteration::getProblem()");
119
120
    if (number == 0) return problem1;
    if (number == 1) return problem2;
Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
    ERROR_EXIT("invalid problem number\n");
    return NULL;
123
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

private:
  ProblemStatBase *problem1;
  ProblemStatBase *problem2;
};
\end{lstlisting}

The class \verb+MyCoupledIteration+ is finished now.

The next class, \verb+Identity+, implements the identity $I(x)=x$ for \verb+double+ variables. An arbitrary degree can be given to it. The class is used later to determine the quadrature degree used for the right hand side of \verb+problem2+. 

\begin{lstlisting}{}
class Identity : public AbstractFunction<double, double>
{
public:
139
  Identity(int degree) : AbstractFunction<double, double>(degree) {}
Thomas Witkowski's avatar
Thomas Witkowski committed
140

141
142
143
144
  double operator()(const double& x) const 
  {
    return x;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
};
\end{lstlisting}

Now, we start with the main program:

\begin{lstlisting}{}
int main(int argc, char* argv[])
{
  FUNCNAME("main");
  TEST_EXIT(argc == 2)("usage: couple initfile\n");
  Parameters::init(true, argv[1]);

  // ===== create and init the first problem ===== 
  ProblemScal problem1("problem1");
  problem1.initialize(INIT_ALL);

  // ===== add boundary conditions for problem1 =====
162
  problem1.addDirichletBC(1, new G);
Thomas Witkowski's avatar
Thomas Witkowski committed
163
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
\end{lstlisting}

So far, we created and initialized \verb+problem1+ and its boundary conditions.

Now, we create \verb+problem2+. It should have its own finite element space, system, solver and file writer, but the mesh should be adopted from \verb+problem1+.  

\begin{lstlisting}{}
  // ===== create and init the second problem ===== 
  Flag initFlag = 
    INIT_FE_SPACE |
    INIT_SYSTEM |
    INIT_SOLVER |
    INIT_FILEWRITER;

  Flag adoptFlag =
    CREATE_MESH |
    INIT_MESH;

  ProblemScal problem2("problem2");
  problem2.initialize(initFlag,
		      &problem1,
		      adoptFlag);
\end{lstlisting}

The operators for the first problem are defined like in Section \ref{s:ellipt code}. Here, we only show the operators of \verb+problem2+.

\begin{lstlisting}{}
  // ===== create operators for problem2 =====
  Operator matrixOperator2(Operator::MATRIX_OPERATOR, 
                           problem2.getFESpace());
193
  matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
196
  problem2.addMatrixOperator(&matrixOperator2);
  
  Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFESpace());
197
198
  rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(),
                                                new Identity(degree)));
Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
201
202
203
204
205
206
207
  problem2.addVectorOperator(&rhsOperator2);
\end{lstlisting}

At the left hand side, we have just an ordinary \verb+Simple_ZOT+. At the right hand side, we have a zero order term of the form $f(u)$ with $f=I$ the identity. $u$ is given by the solution DOF vector of \verb+problem1+. $I$ maps the values of the DOF vector evaluated at quadrature points to itself. The function degree is used to determine the needed quadrature degree in the assembler.
 
Now, the adaptation loop is created:

\begin{lstlisting}{}
  // ===== create adaptation loop and iteration interface =====
208
  AdaptInfo *adaptInfo = new AdaptInfo("couple->adapt", 1);
Thomas Witkowski's avatar
Thomas Witkowski committed
209
210
211

  MyCoupledIteration coupledIteration(&problem1, &problem2);

212
  AdaptStationary *adapt = new AdaptStationary("couple->adapt",
Thomas Witkowski's avatar
Thomas Witkowski committed
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
289
290
291
292
293
294
295
296
                                               &coupledIteration,
                                               adaptInfo);
\end{lstlisting}

Note that not a pointer to one of the problems is passed to the adaptation loop, but a pointer to the \verb+coupledIteration+ object, which in turn knows both problems. 

The adaptation loop is now started. After it is finished, the solutions of both problems are written.

\begin{lstlisting}{}
  // ===== start adaptation loop =====
  adapt->adapt();

  // ===== write solution ===== 
  problem1.writeFiles(adaptInfo, true);
  problem2.writeFiles(adaptInfo, true);
}
\end{lstlisting}

\subsection{Parameter file}

We have one adaptation loop called \verb+couple->adapt+:

\begin{lstlisting}{}
couple->adapt->tolerance:         1e-8
couple->adapt->max iteration:     10
couple->adapt->refine bisections: 2
\end{lstlisting}

The coupled problem consits of two sub problems. The first problem creates the mesh, solves its linear system of equations, estimates the error, adapts the mesh, and finally writes its output:

\begin{lstlisting}{}
coupleMesh->macro file name:    ./macro/macro.stand.2d
coupleMesh->global refinements: 0

problem1->mesh:                 coupleMesh
problem1->dim:                  2
problem1->polynomial degree:    1

problem1->solver:                 cg % no, bicgstab, cg, gmres, odir, ores
problem1->solver->max iteration:  1000
problem1->solver->tolerance:      1.e-8
problem1->solver->left precon:    diag

problem1->estimator:              residual
problem1->estimator->C0:          0.1 % constant of element residual
problem1->estimator->C1:          0.1 % constant of jump residual

problem1->marker->strategy:       2 % 0: no 1: GR 2: MS 3: ES 4:GERS
problem1->marker->MSGamma:        0.5

problem1->output->filename:       output/problem1
problem1->output->AMDiS format:   1
problem1->output->AMDiS mesh ext: .mesh
problem1->output->AMDiS data ext: .dat
\end{lstlisting}

The second problem uses the mesh of \verb+problem1+. So it creates no mesh, no estimator, and no marker. But a solver is needed to solve \verb+problem2+s linear system of equations, and a file writer to write the solution:

\begin{lstlisting}{}
problem2->dim:                  2
problem2->polynomial degree:    1

problem2->solver:                 cg % no, bicgstab, cg, gmres, odir, ores
problem2->solver->max iteration:  1000
problem2->solver->tolerance:      1.e-8
problem2->solver->left precon:    diag

problem2->estimator:              no

problem2->marker->strategy:       0

problem2->output->filename:       output/problem2
problem2->output->AMDiS format:   1
problem2->output->AMDiS mesh ext: .mesh
problem2->output->AMDiS data ext: .dat
\end{lstlisting}

\subsection{Macro file}
We again use the macro file \verb+macro/macro.stand.2d+, which was described in Section \ref{s:ellipt macro}.

\subsection{Output}
The solution of the first problem is written to the files \verb+output/problem1.mesh+ and\\ \verb+output/problem1.dat+. 
The solution of the second problem is written to the files\\ \verb+output/problem2.mesh+ and \verb+output/problem2.dat+. 
We don't visualize the results here, because they conform with the results showed in Section \ref{s:vecellipt output}.