Tutorial/Coupled problems

Revision as of 15:53, 12 December 2011 by Bthierry (talk | contribs) (Mathematics)

Revision as of 15:53, 12 December 2011 by Bthierry (talk | contribs) (Mathematics)

Contents

Introduction

Here is explained how to use the solution of a first problem as a data in a second problem. Two kind of problem are studied here. In the "first" one, the solution of the first problem directly appears in the weak formulation of the second problem (for example, as a source or a Neumann boundary) whereas in the second kind of problem, the solution is used as a constraint (Dirichlet boundary condition), so it should appears in the function space. Note that a Lagrange multiplier could be used to force the Dirichlet constraint in the weak formulation.

First kind of coupled problem : solution appearing in the weak formulation

Mathematics

Let the computation domain be the unit square $\Omega = [0,1]\times[0,1]$ with boundary $\Gamma$ and unit outwardly directed normal $\mathbf{n}$. Consider the two coupled following problems : firstly, find $u$, solution of \begin{equation} \begin{cases}\label{eq:problemU} \Delta u = 0 & \text{in } \Omega,\\ u = C & \text{on }\Gamma, \end{cases} \end{equation} where $C$ is a constant. Then, find the solution $v$ of the second problem \begin{equation} \begin{cases}\label{eq:problemV} -\Delta v + v = 0 & \text{in } \Omega,\\ \displaystyle{\frac{\partial v}{\partial \mathbf{n}} = u} & \text{on }\Gamma, \end{cases} \end{equation} where $u$ is the solution of problem \ref{eq:problemU}. Obviously, the function $u$ is the constant function equal to $C$ and thus, there is no "real" need in solving numerically problem \ref{eq:problemU}, keep in mind that this is just an example. In order to solve numerically these problem using a finite element method, on has to write the weak formulations of problems \ref{eq:problemU} and \ref{eq:problemV}, which read as: \begin{equation}\label{eq:WeakFormulationU} \left\{\begin{array}{l} \text{Find } u\in H^1_C(\Omega) \text{ such that}\\ \displaystyle{\forall u'\in H^1_0(\Omega), \qquad \int_{\Omega} \nabla u\cdot\nabla u' \;{\rm d}\Omega = 0}, \end{array}\right. \end{equation} and \begin{equation}\label{eq:WeakFormulationV} \left\{\begin{array}{l} \text{Find } v\in H^1(\Omega) \text{ such that, }\\ \displaystyle{\forall v'\in H^1(\Omega), \qquad \int_{\Omega} \nabla v\cdot\nabla v' \;{\rm d}\Omega + \int_{\Omega}vv' \;{\rm d}\Omega - \int_{\Gamma}uv'\;{\rm d}\Gamma = 0}, \end{array}\right. \end{equation} where $H^1(\Omega)$ is the classical Sobolev space, $H^1_0(\Omega) = \left\{w\in H^1(\Omega), \text{ such that } w|_{\Gamma}=0\right\}$ and $H^1_C(\Omega) = \left\{w\in H^1(\Omega), \text{ such that } w|_{\Gamma}=C\right\}$.

GetDP

In GetDP, this kind of problem is very easy to solve. Indeed, the user has to introduce two function space (one per solution), and ask GetDP to solve each problem in the right order. There is no trap ! Below are the GMSH and GetDP files. Compared to the other academic examples, the only new things are in the .pro file.

Second kind of coupled problem : solution used as a Dirichlet boundary condition

Mathematics

In this section, the first problem (\ref{eq:problemU}) does not change and $u$ is the solution of this problem, however the second problem (\ref{eq:problemV}) is now to find $v$, solution of \begin{equation} \begin{cases}\label{eq:problemV2} -\Delta v + v = 0 & \text{in } \Omega,\\ \displaystyle{v = u} & \text{on }\Gamma, \end{cases} \end{equation} where $u$ is the solution of problem \ref{eq:problemU}. The solution $u$ is now used as a Dirichlet boundary condition. To obtain the weak formulation of problem (\ref{eq:problemV2}), we introduce the space $H^1_u(\Omega) = \left\{w\in H^1(\Omega) \text{ such that } w|_{\Gamma} = u|_{\Gamma}\right\}$, where $u$ is the solution of problem (\ref{eq:problemU}). Thus, we have \begin{equation}\label{eq:WeakFormulationV2} \left\{\begin{array}{l} \text{Find } v\in H^1_u(\Omega) \text{ such that, }\\ \displaystyle{\forall v'\in H^1_0(\Omega), \qquad \int_{\Omega} \nabla v\cdot\nabla v' \;{\rm d}\Omega + \int_{\Omega}vv' \;{\rm d}\Omega = 0}, \end{array}\right. \end{equation}

GetDP

In GetDP, a Dirichlet boundary condition is generally imposed through the "Constraint" term, it can be imposed with the help of a Lagrange multiplier. In the second case, the Dirichlet boundary condition is no more forced in the function space but appears in the weak formulation. Here is only considered the first case using the GetDP function "AssignFromResolution". The procedure can be summarized as follows:

  • Construct only ONE function space that will contain a data corresponding to $u$ and finally to $v$ ($u$ will not be saved at the end !)
  • Solve problem 1 (and save $u$)
  • Transfert solution of problem 1 to problem 2
  • Solve problem 2 (and erase $u$ by $v$)

Sometimes it is usefull to save both $u$ and $v$ at the end of the calculation. This can be done by adding an auxiliary problem that "copy" $u$ (see the last paragraph).

Files

Saving the solution $u$