couple.tex 12.1 KB
 Thomas Witkowski committed Oct 23, 2008 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 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 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 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 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 289 290 291 292 293 294 295 296 297 298 299 \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 v = u. 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) {}; \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"); }; void endIteration(AdaptInfo *adaptInfo) { FUNCNAME("StandardProblemIteration::endIteration()"); MSG("\n"); MSG("end of iteration %d\n", adaptInfo->getSpaceIteration()+1); MSG("=============================\n"); }; \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; if(toDo.isSet(MARK)) markFlag = problem1->markElements(adaptInfo); if(toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED)) { flag = problem1->refineMesh(adaptInfo); } if(toDo.isSet(BUILD)) problem1->buildAfterCoarsen(adaptInfo, markFlag); if(toDo.isSet(SOLVE)) problem1->solve(adaptInfo); if(toDo.isSet(BUILD)) problem2->buildAfterCoarsen(adaptInfo, markFlag); if(toDo.isSet(SOLVE)) problem2->solve(adaptInfo); if(toDo.isSet(ESTIMATE)) problem1->estimate(adaptInfo); return flag; }; \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; }; ProblemStatBase *getProblem(int number = 0) { FUNCNAME("CoupledIteration::getProblem()"); if(number == 0) return problem1; if(number == 1) return problem2; ERROR_EXIT("invalid problem number\n"); return NULL; }; 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 { public: MEMORY_MANAGED(Identity); Identity(int degree) : AbstractFunction(degree) {}; const double& operator()(const double& x) const { static double result; result = x; return result; }; }; \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 ===== problem1.addDirichletBC(1, NEW G); \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()); matrixOperator2.addZeroOrderTerm(NEW Simple_ZOT); problem2.addMatrixOperator(&matrixOperator2); Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFESpace()); rhsOperator2.addZeroOrderTerm(NEW VecAtQP_ZOT(problem1.getSolution(), NEW Identity(degree))); 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 ===== AdaptInfo *adaptInfo = NEW AdaptInfo("couple->adapt", 1); MyCoupledIteration coupledIteration(&problem1, &problem2); AdaptStationary *adapt = NEW AdaptStationary("couple->adapt", &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}.