text:
\subsection{Discretized coupling formulation and coupled problem}
LP{} conservation equations \cref{eq:0D_mass_equation,eq:0D_temperature_equation} are then discretized in time at a lower level with a method of choice (e.g.\ explicit, implicit, Euler, Runge-Kutta, multistep methods, etc.) and then coupled in time on a higher level resulting in a two-level time scheme. At the highest level, those discretized equations are solved and coupled between times $t^0,t^1,\ \ldots,t^n,\ \ldots$ with a macro time step~$\Delta t$ and synchronized at each of these times. At the lowest level, each subdomain manages its own time integration scheme, its own micro time step $\delta t$ and its own internal time. The integration scheme used by the subdomain is assumed to be adapted to the physical local problem. In the following, discretized values evaluated at time $t^n$ are denoted by the superscript $n$.

We adopt the following notations. For each coupled subdomain $\{ \domain{i} \}_{i\in[1,m]}$, discretized equations are represented by a function $\mathcal{F}_i^{\Delta t}$ used to solve and advance in time the problem of one time step $\Delta t$. This function takes as parameters the state vector $\mathbf{u}_i$ and the input interface variables $\{\mathbf{b}_{ji}\}_{j \in N_i} \overset{\text{def.}}{=} \mathbf{b}_{\star i}$ of domain $\domain{i}$. Equations of the coupled problem are then given by
\begin{equation}
\label{eq:flow_formulation}
\left\{
	\begin{aligned}
		\mathcal{F}_1^{\Delta t}(\mathbf{u_1},&\ \{\mathbf{b}_{j1}(\mathbf{u}_j,\ \mathbf{b}_{\star j})\}_{j \in N_1}) = 0\\
		\mathcal{F}_2^{\Delta t}(\mathbf{u_2},&\ \{\mathbf{b}_{j2}(\mathbf{u}_j,\ \mathbf{b}_{\star j})\}_{j \in N_2}) = 0\\
		&\vdots\\
		\mathcal{F}_m^{\Delta t}(\mathbf{u_m},&\ \{\mathbf{b}_{jm}(\mathbf{u}_j,\ \mathbf{b}_{\star j})\}_{j \in N_m}) = 0\\
	\end{aligned}
\right.
\end{equation}
in such a way that the inter dependencies between domain $\domain{i}$ and its neighbors are highlighted by the closure law functions $\{\mathbf{b}_{ji}\}_{j \in N_i}$ which is function of the coupled subdomain internal variable $\mathbf{u}_j$ and interface variables $\mathbf{b}_{\star j}$, which might in turn need values from domain $\domain{i}$ to be calculated.

Those interface dependencies highlight the input-output relationship between interface variables of the different subdomains. This suggests that solvers can be seen as closed entities or ``black-boxes'' taking input interface variables $\mathbf{b}_{\star i}$ from neighboring subdomains and giving back output interface variables $\mathbf{b}_{i\star}$ to neighboring subdomains. Therefore, it seems natural to represent solver of domain $\domain{i}$ as a function $\mathcal{M}_i^{\Delta t} : \mathbb{R}^{3 \times n_i} \mapsto \mathbb{R}^{3 \times n_i}$ in which we explicitly hide subdomain state variable $\mathbf{u}_i$ and its integration, and only the input and output interface variables are visible. It is possible that we have access to limited information about them, e.g.\ no derivative data. For instance, solver of domain $\domain{i}$ is defined by
\begin{equation}
\label{eq:solver_formulation}
\mathbf{b}_{i\star} = \mathcal{M}_i^{\Delta t}(\mathbf{b}_{\star i}).
\end{equation}

From \cref{eq:solver_formulation}, we can define the coupled problem at interface $\interface{ij}$ in term of solvers :
\begin{align}
\mathbf{b}_{ij} &= \mathcal{P}_{ij} \circ \mathcal{M}_i^{\Delta t}\left(\mathbf{b}_{ji},\{\mathbf{b}_{ki}\}_{k \in N_i,\ k \neq j}\right)\label{eq:ij_projection},\\
\mathbf{b}_{ji} &= \mathcal{P}_{ji} \circ \mathcal{M}_j^{\Delta t}\left(\mathbf{b}_{ij},\{\mathbf{b}_{kj}\}_{k \in N_j,\ k \neq i}\right)\label{eq:ji_projection}
\end{align}
with $\mathcal{P}_{ij}$ the interface $\interface{ij}$ projector. Those equations emphasize the ``action-reaction'' at interface between the two domains : a slight modification of interface variable $\mathbf{b}_{ji}$ will in turn change the associated interface variable $\mathbf{b}_{ij}$, and vice versa. The strength of the associated coupling is represented by the Jacobian matrices $\{\frac{\partial \mathcal{M}_i^{\Delta t}}{\partial \mathbf{b}_{ki}}\}_{k \in N_i}$ and $\{\frac{\partial \mathcal{M}_j^{\Delta t}}{\partial \mathbf{b}_{kj}}\}_{k \in N_j}$. Moreover it is difficult to evaluate this strength since they may not be available.

Combining \cref{eq:ij_projection} with \cref{eq:ji_projection} gives the widely used fixed point equation at interface $\interface{ij}$ (see \cite{mehl_parallel_2016})
\begin{equation}
\label{eq:fpe_formulation}
\mathbf{b}_{ij} = \mathcal{P}_{ij} \circ \mathcal{M}_i^{\Delta t}\left(\mathcal{P}_{ji} \circ \mathcal{M}_j^{\Delta t}\left(\mathbf{b}_{ij},\{\mathbf{b}_{kj}\}_{k \in N_j,\ k \neq i}\right),\{\mathbf{b}_{ki}\}_{k \in N_i,\ k \neq j}\right)
\end{equation}
which allows us to define the residual operator $\mathcal{R}_{ij}$ at interface $\interface{ij}$  by
\begin{equation}
\label{eq:residual_formulation}
\mathcal{R}_{ij}(\mathbf{b}_{ij}) \overset{\textit{def}}{=} \mathcal{P}_{ij} \circ \mathcal{M}_i^{\Delta t}\left(\mathcal{P}_{ji} \circ \mathcal{M}_j^{\Delta t}\left(\mathbf{b}_{ij},\{\mathbf{b}_{kj}\}_{k \in N_j,\ k \neq i}\right),\{\mathbf{b}_{ki}\}_{k \in N_i,\ k \neq j}\right) - \mathbf{b}_{ij}
\end{equation}
for a guess candidate $\mathbf{b}_{ij}$. The residual $\mathcal{R}_{ij}(\mathbf{b}_{ij})$ will express the imbalance created at interface $\interface{ij}$ : if domains $\domain{i}$ and $\domain{j}$ are strongly coupled, equilibrium conditions at interface are fulfilled and the interface residual is null, otherwise the two domains are only weakly coupled.

\subsection{Coupling schemes}
\label{sec:coupling_scheme}
\paragraph{Explicit coupling schemes (ECS{}s).}Referred to as ``conventional serial staggered'' in \cite{felippa_staggered_1980}, ECS{}s solve the coupled problem with one call of solver per time step. They are based on a Gauss-Seidel semi-explicit solution of the coupled problem \cref{eq:flow_formulation}. Obviously, another scheme based on a fully-explicit or Jacobi solution can be used to allow for more algorithm parallelism. If we assume that solver $\mathcal{M}_i^{\Delta t}$ is solved before solver $\mathcal{M}_j^{\Delta t}$, the scheme solves at interface $\interface{ij}$ :
\begin{align}
\mathbf{b}_{ij}^{n+1} &= \mathcal{P}_{ij} \circ \mathcal{M}_i^{\Delta t}\left(\mathbf{b}_{ji}^{n},\{\mathbf{b}_{ki}^{\bullet}\}_{k \in N_i,\ k \neq j}\right),\\
\mathbf{b}_{ji}^{n+1} &= \mathcal{P}_{ji} \circ \mathcal{M}_j^{\Delta t}\left(\mathbf{b}_{ij}^{n+1},\{\mathbf{b}_{kj}^{\bullet}\}_{k \in N_j,\ k \neq i}\right)
\end{align}
with $\mathbf{b}_{ki}^{\bullet}$ is evaluated at $t^n$ or $t^{n+1}$ if solver $\mathcal{M}_k^{\Delta t}$ is called before or after solver $\mathcal{M}_i^{\Delta t}$.

While potentially attractive and fast since only one call per solver is done during each time step, it is well known that this scheme yields poor accuracy and stability issues \cite{piperno_partitioned_1995}. Because of the \textit{time-lag} caused by the semi-explicit resolution, it is unlikely that interface residuals defined by \cref{eq:residual_formulation} are null and equilibrium conditions are not enforced at interface. Besides, despite several improvements and studies in \cite{farhat_mixed_1995} \cite{guillard_significance_2000} \cite{farhat_consistency_1990} \cite{farhat_robust_2010} \cite{farhat_unconditionally_1991} \cite{piperno_partitioned_2001}, the weak coupling reached by explicit coupling schemes is often not enough and only a strong coupling at the end of the time step can ensure proper stability properties \cite{causin_added-mass_2005} \cite{giles_stability_1997}.

\paragraph{Implicit coupling schemes (ICS{}s).} A way to fix the \textit{time-lag} is to use implicit coupling between subdomains. At interface $\interface{ij}$, implicit coupling gives~:
\begin{align}
\mathbf{b}_{ij}^{n+1} &= \mathcal{P}_{ij} \circ \mathcal{M}_i^{\Delta t}\left(\mathbf{b}_{ji}^{n+1},\{\mathbf{b}_{ki}^{n+1}\}_{k \in N_i,\ k \neq j}\right)\label{eq:implicit_coupling_1},\\
\mathbf{b}_{ji}^{n+1} &= \mathcal{P}_{ji} \circ \mathcal{M}_j^{\Delta t}\left(\mathbf{b}_{ij}^{n+1},\{\mathbf{b}_{kj}^{n+1}\}_{k \in N_j,\ k \neq i}\right)\label{eq:implicit_coupling_2}
\end{align}

We clearly see that \cref{eq:implicit_coupling_1,eq:implicit_coupling_2} give discrete interface values respecting the fixed point \cref{eq:fpe_formulation} leading to a null interface residual. Hence, interface $\interface{ij}$ is at equilibrium and a strong coupling between equations is obtained. However, while being mathematically interesting, \cref{eq:implicit_coupling_1,eq:implicit_coupling_2} do not allow to decouple the two equations at each time step and it is often more convenient to use iterative methods to solve them. By iterative methods we mean methods using iterative processes inside one coupling iteration, i.e. between times $t^n$ and $t^{n+1}$, until a convergence criterion is reached. In that case interface variables \textit{verify} (with an error bounded by the convergence criterion) at each interface the associated fixed point \cref{eq:fpe_formulation} leading to a strong coupling between equations of the coupled problem. At interface $\interface{ij}$ between domain $\domain{i}$ and $\domain{j}$, the block Gauss-Seidel coupling scheme, also known as staggered or partitioned coupling scheme, is an iterative method based on Newton-Raphson iterations of the fixed point \cref{eq:fpe_formulation} given by
\begin{equation}
\label{eq:newton_iterations}
\mathbf{b}_{ji}^{n+1, k+1}  = \mathbf{b}_{ji}^{n+1, k}  + \left(\frac{\mathrm{d}\mathcal{R}_{ji}}{\mathrm{d}\mathbf{b}_{ji}}\Bigr|_{\mathbf{b}_{ji}^{n+1, k}}\right)^{-1} (-\mathcal{R}_{ji}(\mathbf{b}_{ji}^{n+1,k}))
\end{equation}
with the residual $\mathcal{R}_{ji}(\mathbf{b}_{ji}^{n+1, k})$ at iteration $k$ defined by
\begin{equation}
\label{eq:residual}
\mathcal{R}_{ji}(\mathbf{b}_{ji}^{n+1, k}) = \mathbf{r}_{ji}^{n+1,k} = \mathcal{M}_j^{\Delta t} \circ \mathcal{M}_i^{\Delta t}(\mathbf{b}_{ji}^{n+1, k}) - \mathbf{b}_{ji}^{n+1, k}
\end{equation}
which allows us to define a convergence criterion for the iterative process by
\begin{equation}
\label{eq:convergence}
\frac{\norm{\mathcal{M}_j^{\Delta t} \circ \mathcal{M}_i^{\Delta t}(\mathbf{b}_{ji}^{n+1, k}) - \mathbf{b}_{ji}^{n+1, k}}}{\norm{\mathbf{b}_{ji}^{n+1, k}}} = \frac{\norm{\mathbf{r}_{ji}^{n+1,k}}}{\norm{\mathbf{b}_{ji}^{n+1, k}}} \leq \epsilon_{\textrm{rel}}
\end{equation}
with $\epsilon_{\textrm{rel}}$ the relative stopping criterion for the subiterative process. \textit{Gauss-Seidel} or \textit{staggered} symbolizes the way interfaces data are sequenced through solvers (or \textit{block}) inside coupling iterations.
\paragraph{Remark on Jacobian matrices and associated solvers.} In \cref{eq:newton_iterations}, the Jacobian matrix or its inverse are usually not known and/or not calculable since solvers are seen as black-box solvers, thus the iterations can only be approximated by approximation of the Jacobian matrix or its inverse leading to Quasi-Newton techniques. In \cite{gerbeau_quasi-newton_2003}, Gerbeau and al. use reduced order models to calculate the Jacobian matrix. In \cite{michler_interface_2005}, Michlet and al. use Newton-Krylov subiterations to approximate the Jacobian matrix leading to the so-called \textit{Jacobian free Newton-Krylov} method. From previous residuals, Degroote and al. in \cite{degroote_performance_2009} solve least-squares problems to approximate the inverse of the Jacobian matrix leading to \textit{Interface Quasi-Newton - Inverse Least-Squares} (IQN-ILS) methods while Vierendeels and al. \cite{vierendeels_implicit_2007} approximate the Jacobian matrix with similar methods leading to \textit{Interface Quasi-Newton - Least-Squares} (IQN-LS). Multigrid Quasi-Newton methods are also studied in \cite{degroote_multi-level_2012}. A recent review of different Quasi-Newton methods can be found in \cite{minami_performance_2010} for fluid-structure interaction and in \cite{ganine_nonlinear_2013} for thermal-structure interaction.

Much cheaper methods use relaxation techniques for the interface variable iterations which consists in using an approximation of the inverse of the interface residual Jacobian matrix of the form
\begin{equation}
\left(\frac{\mathrm{d}\mathcal{R}_{ji}}{\mathrm{d}\mathbf{b}_{ji}}\Bigr|_{\mathbf{b}_{ji}^{n+1, k}}\right)^{-1} \approx -\omega^k\,I
\end{equation}
for which the fixed point iterations with dynamic relaxation are given by
\begin{equation}
\label{eq:relaxation}
\mathbf{b}_{ji}^{n+1,k+1} = \omega^k \mathcal{M}_j^{\Delta t} \circ \mathcal{M}_i^{\Delta t}(\mathbf{b}_{ji}^{n+1, k}) + (1-\omega^k) \mathbf{b}_{ji}^{n+1, k}.
\end{equation}
Note that for $\omega^k=1$, iterations given by \cref{eq:relaxation} are classical Picard iterations which converge generally only linearly and very slowly. In \cite{kuttler_fixed-point_2008}, the authors show that the block Gauss-Seidel iterative method with relaxation techniques for the resolution of equation \cref{eq:newton_iterations} can be very efficient at a surprisingly low cost in comparison to more elaborated Quasi-Newton methods. Besides, it is shown in \cite{ramiere_iterative_2015} that these techniques are also proven to be very competitive in comparison to the direct solution of the nonlinear problem given by \cref{eq:newton_iterations} with traditional Newton-Raphson methods using derivative data. In this paper, we focus on block Gauss-Seidel iterative method with relaxation techniques.

For instance, the solution of equation \cref{eq:newton_iterations} by the Steffensen's method accelerated with Aitken's delta-squared $\Delta^2$ method \cite{aitken_xxv.bernoullis_1927} can be cast into fixed point iterations with dynamic relaxation leading to methods of order $2$, see \cite{ramiere_iterative_2015}. Less costly but with a convergence rate of the golden ratio $\frac{1+\sqrt{5}}{2}$, the fixed point iterations given by the secant method can also be seen as fixed point iterations with dynamic relaxation leading to
\begin{equation}
\omega^k = -\omega^{k-1}\frac{\inner{\mathbf{r}_{ji}^{n+1,k}-\mathbf{r}_{ji}^{n+1,k-1}}{\mathbf{r}_{ji}^{n+1,k-1}}}{\inner{\mathbf{r}_{ji}^{n+1,k}-\mathbf{r}_{ji}^{n+1,k-1}}{\mathbf{r}_{ji}^{n+1,k}-\mathbf{r}_{ji}^{n+1,k-1}}}.
\end{equation}
One can also use constant relaxation given by a constant parameter $\omega$. However, this requires the determination of the best relaxation parameters, i.e. leading to the highest rates of convergence, which is highly problem-dependent.

\paragraph{Implicit coupling schemes for event detection and model synchronization.} In LP{} modeling, state transitions are triggered by internal events corresponding to activation of threshold functions. This activation is dependent on input from the coupled models and thus times of events are unknowns of the coupled problem. As depicted \cref{fig:domain_states}, each state has its own solver, its own set of equations and its own interface and boundary conditions.
\begin{figure}
    \centering
	\includegraphics[width=0.7\textwidth]{domain_states.eps}
	\caption{The abstract representation of states $\domain[(1)]{i}$ and $\domain[(2)]{i}$ of domain $\domain{i}$. Each state has its own interfaces $\interface[(1)]{\bullet}$ and $\interface[(2)]{\bullet}$ with associated variables $\mathbf{b}^{(1)}_{\bullet}$ and $\mathbf{b}^{(2)}_{\bullet}$. A state transition occurs when a transition function crosses zero.}
	\label{fig:domain_states}
\end{figure}

State transition can lead to discontinuities of state or interface variables. A missed event can lead to inconsistent and incoherent physical states, e.g. the disappearance of a model not seen by other coupled models, and/or mathematical difficulties with non-smooth functions in the coupled problem \cref{eq:flow_formulation} bringing usual theorems out of their scope of validity, with possibly breakdown issue like divergence of Newton iterations.

During one coupling time step $\Delta t$, the coupling scheme must be able to adapt its time step to synchronize all the models on the time of the first event. However, synchronization strategies are hard to achieve with ECS{}s which can often only synchronized coupled models at each coupled time step leading to time detection errors of order $\mathcal{O}(\Delta t)$. Hereafter, we describe how the previous ICS{}s can be used for event detection to ensure proper synchronization between models.\medskip

Let us consider the coupling at interface $\interface{ij}$ between domain $i$ and $j$. As depicted in \cref{fig:events_time_loop}, internal events are triggered during one coupling time step $\Delta t$ by solver $\mathcal{M}_i$ and/or solver $\mathcal{M}_j$ at respective times $\textcolor{red!50!black}{t_i^{\star}}$ and $\textcolor{red!50!black}{t_j^{\star}}$ between $t^n$ and $t^n+\Delta t$. When an event occurs, the solver stops its computation and do not compute the whole coupling time step $\Delta t$. Solvers should be synchronized on the first triggered event at time $\textcolor{red!50!black}{t^{\star}} = \textcolor{red!50!black}{t_j^{\star}}$.
\begin{figure}
\begin{tikzpicture}
\node at (-2.3,0.) {$\mathcal{M}_i$ \textbf{internal} time loop};
\draw [->, thick] (0,0) --  (8.3,0);
\node at (-2.3,-1.0) {$\mathcal{M}_j$  \textbf{internal} time loop};
\draw [->, thick] (0,-1) --  (8.3,-1);
\node at (-2.3,-2.0) {\textbf{Shared} time loop};
\draw [->, thick] (0,-2) --  (8.3,-2);
\draw[thick] (0.1,-0.1) -- (0.1,0.1);
\draw[thick] (0.1,-1.1) -- (0.1,-0.9);
\draw[thick] (0.1,-2.1) -- (0.1,-1.9);
\node at (0.1,-2.5) {$t^n$};
\draw[thick] (8.0,-0.1) -- (8.0,0.1);
\draw[thick] (8.0,-1.1) -- (8.0,-0.9);
\draw[thick] (8.0,-2.1) -- (8.0,-1.9);
\node at (8.2,-2.5) {$t^n+\Delta t$};
\node at (6.8,-1.6) {$\textcolor{blue!50!black}{t^{\star,k}}$};
\draw (6.7,-1.9) -- (6.9,-2.1);
\draw (6.7,-2.1) -- (6.9,-1.9);
\node at (6,0.4) {$\textcolor{red!50!black}{t_i^{\star}}$};
\draw (5.9,0.1) -- (6.1,-0.1);
\draw (5.9,-0.1) -- (6.1,0.1);
\node at (5.2,-0.6) {$\textcolor{red!50!black}{t_j^{\star}}$};
\draw (5.1,-0.9) -- (5.3,-1.1);
\draw (5.1,-1.1) -- (5.3,-0.9);
\node at (5.2,-1.6) {$\textcolor{red!50!black}{t^{\star}}$};
\draw (5.1,-1.9) -- (5.3,-2.1);
\draw (5.1,-2.1) -- (5.3,-1.9);
\node at (3.0,0.4) {$\textcolor{blue!50!black}{t_i^{\star,k}}$};
\draw (2.9,0.1) -- (3.1,-0.1);
\draw (2.9,-0.1) -- (3.1,0.1);
\node at (4.2,-0.6) {$\textcolor{blue!50!black}{t_j^{\star,k}}$};
\draw (4.1,-0.9) -- (4.3,-1.1);
\draw (4.1,-1.1) -- (4.3,-0.9);
\node at (2.0,-1.6) {$\textcolor{blue!50!black}{\tilde{t}^{\star,k}}$};
\draw (1.9,-1.9) -- (2.1,-2.1);
\draw (1.9,-2.1) -- (2.1,-1.9);
\node at (3.5,-1.6) {$\textcolor{blue!50!black}{t^{\star,k+1}}$};
\draw (3.4,-1.9) -- (3.6,-2.1);
\draw (3.4,-2.1) -- (3.6,-1.9);
\node at (3.5,-2.4) {$\scriptstyle\alpha \textcolor{blue!50!black}{\tilde{t}^{\star,k}} + (1-\alpha)\textcolor{blue!50!black}{t^{\star,k}}$};
\end{tikzpicture}
\caption{Shared coupling time loop and internal time loops of solvers $\mathcal{M}_i$ and $\mathcal{M}_j$. Internal events at times $\textcolor{red!50!black}{t_i^{\star}}$ and $\textcolor{red!50!black}{t_j^{\star}}$ are to be detected by the scheme in order to synchronize on time $\textcolor{red!50!black}{t^{\star}}$ of the first event. An ICS{} detects event at times $\textcolor{blue!50!black}{t_i^{\star,k}}$ and $\textcolor{blue!50!black}{t_j^{\star,k}}$ and has to calculate the new iteration $\textcolor{blue!50!black}{t^{\star,k+1}}$ in such a way that $\textcolor{blue!50!black}{t^{\star,\infty}} \approx \textcolor{red!50!black}{t^{\star}}$.}
\label{fig:events_time_loop}
\end{figure}

During a fixed point iteration of an ICS{} starting at time $t^n$ and ending at the previously calculated time of first event $\textcolor{blue!50!black}{t^{\star,k}}$, events are trigerred at times $\textcolor{blue!50!black}{t_i^{\star,k}}$ and $\textcolor{blue!50!black}{t_j^{\star,k}}$ by solvers. The aim is to build an iterative algorithm calculating the new iteration $\textcolor{blue!50!black}{t^{\star,k+1}}$ in such a way that~: $\textcolor{blue!50!black}{t^{\star,\infty}} \approx \textcolor{red!50!black}{t^{\star}}$. \Cref{alg:event_time_iteration} gives a sequence to calculate this new iteration.
\begin{algorithm}
\caption{Calculate the new iteration $\textcolor{blue!50!black}{t^{\star,k+1}}$}
\label{alg:event_time_iteration}
\renewcommand{\algorithmiccomment}[1]{/* #1 */}
\algsetup{indent=2em}
\begin{algorithmic}[1]
\STATE Compute a new fixed point iteration between solvers $\mathcal{M}_i^{\Delta t}$ and $\mathcal{M}_j^{\Delta t}$ with an ICS{};
\IF {\textit{an event has occurred for model $l\in\{i,j\}$}}
        \STATE $\textcolor{blue!50!black}{t_l^{\star,k}}\gets$ the corresponding event time;
\ELSE
        \STATE $\textcolor{blue!50!black}{t_l^{\star,k}} \gets t^n+\Delta t$;
		\STATE \COMMENT{the whole coupling time step $\Delta t$ was computed.}
\ENDIF
\STATE $\textcolor{blue!50!black}{\tilde{t}^{\star,k}} \gets \min{}(\textcolor{blue!50!black}{t_i^{\star,k}},\textcolor{blue!50!black}{t_j^{\star,k}});$
\IF[$\epsilon_{\textrm{rel}}$ being a tolerance given by the ICS{}.] {$(\vert\textcolor{blue!50!black}{\tilde{t}^{\star,k}} - \textcolor{blue!50!black}{t^{\star,k}}\vert / \Delta t < \epsilon_{\textrm{rel}}$) \AND (\textit{the ICS{} has converged at interface} $\interface{ij}$)}
	\STATE $\textcolor{blue!50!black}{t^{\star,k+1}} \gets \textcolor{blue!50!black}{\tilde{t}^{\star,k}};$
	\STATE $\textcolor{blue!50!black}{t^{\star,\infty}} \gets \textcolor{blue!50!black}{t^{\star,k+1}};$
\ELSE
	\STATE $\textcolor{blue!50!black}{t^{\star,k+1}} \gets \alpha \textcolor{blue!50!black}{\tilde{t}^{\star,k}} + (1-\alpha)\textcolor{blue!50!black}{t^{\star,k}};$
	\STATE \COMMENT{$\alpha\in]0,1[$ being a relaxation parameter given by the ICS{}.}
\ENDIF
\end{algorithmic}
\end{algorithm}
At convergence, the coupled problem has been solved between times $t^n$ and $t^{n+1} = \textcolor{blue!50!black}{t^{\star,\infty}}$. The next coupling time step then starts at time $t^{n+1}$ with the same algorithm.
\paragraph{Concluding remarks regarding the designed iterative algorithm.}At convergence $k=\infty$ of a coupling time step between times $t^n$ and $t^{n+1}$ at interface $\interface{ij}$, we get that :
\begin{itemize}
\item Model $i$ and model $j$ are strongly coupled because the fixed point \cref{eq:fpe_formulation} is verified with at most an error bounded by the tolerance $\epsilon_{\textrm{rel}}$ of the scheme, i.e. \[\frac{\norm{\mathcal{M}_j^{\Delta t} \circ \mathcal{M}_i^{\Delta t}(\mathbf{b}_{ji}^{n+1, \infty}) - \mathbf{b}_{ji}^{n+1, \infty}}}{\norm{\mathbf{b}_{ji}^{n+1, \infty}}} \leq \epsilon_{\textrm{rel}}.\]
\item Model $i$ and model $j$ are synchronized on potential internal triggered events. In particular, they are synchronized on the time of the first event with an error also bounded by the tolerance $\epsilon_{\textrm{rel}}$, i.e. \[\frac{\abs{t^{n+1} - t^{\star}}}{\Delta t} = \frac{\abs{t^{\star, \infty} - t^{\star}}}{\Delta t} \leq  \epsilon_{\textrm{rel}}\]
\end{itemize}

\subsection{Subdomain lumped parameter equations}
\label{sec:subdomain_equations}
Subdomain $\domain{i}$ equations are expressed in terms of mass and energy macroscopic conservation equations (momentum conservation equations are modeled by closure laws, e.g. correlations for heat transfer coefficient). They are obtained from the local conservation equations, here Navier-Stokes equations under the Boussinesq approximation for liquid domains and heat equations for solid domains, integrated over the corresponding subdomain (see \cite{le_tellier_treatment_2017}). Tightly linked to the local physical model, this approach leads to the so-called LP{} model or ``0D'' model of the subdomain  described by the two ordinary differential equations (ODE)
\begin{align}
\D{\mass{i}}{t}  &= \sum_{j \in N_i} \mFlowRate{ij}\quad\text{in}\;\domain{i} \label{eq:0D_mass_equation},\\
\mass{i} \heatCp{i} \D{\temperature{i}}{t} + \sum_{j \in N_i} \mFlowRate{ij} \heatCp{i} (\temperature{i} - \temperature{ij})&= \sigma_i \heatFlux{i} \area{i} + \sum_{j \in N_i} \sigma_{ij} \heatFlux{ij} \area{ij} + \mass{i} \massPower{i}\quad\text{in}\;\domain{i} \label{eq:0D_temperature_equation}
\end{align}
with $\mass{i}\,[\text{kg}{}]$ the mass and $\temperature{i}\,[\text{K}{}]$ the average temperature of $\domain{i}$, $\heatFlux{i}\,[\text{W.m$^{-2}$}{}]$ the heating ($\sigma_{i} = 1$) or cooling ($\sigma_{i} = -1$) heat flux through boundary $\interface{i} = \frontier{} \cap \frontier{i}$ with temperature $\temperature{b_i}$ and area $\area{i}\,[\text{m$^2$}{}]$, $\heatFlux{ij}$ the heating or cooling heat flux and $\mFlowRate{ij}\,[\text{kg.s$ ^{-1}$}{}]$ the algebraic mass flow rate through $\interface{ij}$ with temperature $\temperature{ij}$ and area $\area{ij}$. Finally, $\heatCp{i}\,[\text{J.kg$^{-1}$.K$^{-1}$}{}]$ is the heat capacity and $\massPower{i}\,[\text{W.kg$^{-1}$}{}]$ is the residual power per mass unit coming from fission products of subdomain $\domain{i}$.

The previous physical parameters are obtained from closure laws described hereafter in \cref{sec:interface_equations}. In particular, geometry dependent values, i.e.\ the characteristic length of a domain or the surface of an interface or the volume of a domain, are given by algebraic geometry equations. They take the form of algebraic functions, e.g.\ for the previous surface :
\begin{equation}
\label{eq:0D_geometry}
\area{ij} = \area{ij}(\massDensity{i},\length{i},\vol{i})
\end{equation}
with $\length{i}\,[\text{m}{}]$ the characteristic length and $\vol{i}\,[\text{m$^3$}{}]$ the volume of domain $\domain{i}$. Altogether, ordinary differential \cref{eq:0D_mass_equation,eq:0D_temperature_equation} combined with expressions like \cref{eq:0D_geometry} can be considered as \textit{Differential Algebraic Equations} (DAE's) \cite{petzold_differential/algebraic_1982} describing the LP{} model.

\subsection{Interface equations}
\label{sec:interface_equations}
In the present article, interfaces between two subdomains $\domain{i}$ and $\domain{j}$ can be of two types : free moving or fixed boundaries. Physically moving boundaries are, for example, boundaries between a liquid and a solid domain exchanging melted or solid materials. They are associated with a plane fusion solidification front corresponding to the Stefan condition at the interface (see \cite{le_tellier_treatment_2017} for further details). Fixed interfaces correspond to thermal equilibrium assuming no mass exchange through the interface and thermal conduction.

In case of a mobile interface, equilibrium conditions at the interface $\interface{ij}$ are given by
\begin{align}
\label{eq:fusion_front_equilibrium_condition_1}
\mFlowRate{ij} &= -\mFlowRate{ji}\quad\text{in}\;\interface{ij},\\
\label{eq:fusion_front_equilibrium_condition_2}
\heatFlux{ij}\area{ij} &= - \heatFlux{ji}\area{ji} + \dEnthalpy[\textrm{fus.}]{} \mFlowRate{ij}\quad\text{in}\;\interface{ij},\\
\label{eq:fusion_front_equilibrium_condition_3}
\temperature{ij} &= \temperature{ji} = \temperature[\textrm{fus.}]{}\quad\text{in}\;\interface{ij}
\end{align}
with $\dEnthalpy[\textrm{fus.}]{}$ the fusion enthalpy and $\temperature[\textrm{fus.}]{}$ the fusion temperature of domain $\domain{i}$, both assumed to be fixed. In particular, those conditions stipulate that the mass flow rate should be the same on both sides of the interface for mass conservation and that the heat fluxes should respect the Stefan condition. To simplify subdomain materials are treated as pure body and no thermochemistry is considered.

In case of a fixed interface, the thermal equilibrium conditions at interface $\interface{ij}$ are given by
\begin{align}
\label{eq:thermal_equilibrium_condition_1}
\mFlowRate{ij} &= -\mFlowRate{ji} = 0\quad\text{in}\;\interface{ij},\\
\label{eq:thermal_equilibrium_condition_2}
\temperature{ij} &= \temperature{ji}\quad\text{in}\;\interface{ij},\\
\label{eq:thermal_equilibrium_condition_3}
\heatFlux{ij}\area{ij} &= -\heatFlux{ji}\area{ji}\quad\text{in}\;\interface{ij}
\end{align}

Obviously, with the previous equations, appropriate closure laws for interface heat fluxes and temperatures are required. In traditional mesh based models, interface variables are given by a projection operator which are coherent with the domain equations, e.g.\ the restriction of the variables values over all the domain $\domain{i}$ to its interface $\interface{ij} = \frontier{i} \cap \frontier{j}$, and thus are consistent with the physical equations. For LP{} models, interface variable $\mathbf{b}_{ij}$ is calculated from spatially averaged data $\mathbf{u}_i$ from domain $\domain{i}$ and data $(\{\mathbf{b}_{ij}\}_{j \in N_i})  \overset{\textit{def}}{=} \mathbf{b}_{i\star}$ from the interfaces. They are given by closure laws which take the form of algebraic functions $\mathbf{b}_{ij} = \mathbf{b}_{ij}(\mathbf{u}_i, \mathbf{b}_{i\star})$. Thus they propagate data instantly, from one specific interface to all the other ones.

If the domain $\domain{i}$ is solid, such closure law functions for heat fluxes $\heatFlux{i\bullet}$ can be calculated from the heat diffusion conduction equation under certain assumptions and approximations. A comparison of those different approximate models with a reference solution given by a finite element discretization of the heat conduction equation can be found in \cite{le_tellier_treatment_2017}. For example, the \textit{stationary} model uses a 1D cylindrical with adiabatic lateral boundaries approximation and assumes a quadratic temperature profile in the solid domain. Conduction heat fluxes at interfaces $\interface{ij}$ and $\interface{ik}$ associated with the upper and lower cylinder surfaces are then given by the temperature derivative w.r.t the spatial direction leading to the closure law functions \cite{le_tellier_treatment_2017} :
\begin{align}
\heatFlux{ij} &= \heatFlux{ij}(\mathbf{u}_i, \mathbf{b}_{ik}) = \thermalCond{i} \frac{6 \temperature{i} - 4 \temperature{ij} - 2 \temperature{ik}}{\length{i}}, \label{eq:heat_flux_stationary_model_1}\\
\heatFlux{ik} &= \heatFlux{ik}(\mathbf{u}_i, \mathbf{b}_{ij}) = \thermalCond{i} \frac{6 \temperature{i} - 4 \temperature{ik} - 2 \temperature{ij}}{\length{i}}. \label{eq:heat_flux_stationary_model_2}
\end{align}
Note, for instance, the instant propagation of interface $\interface{ik}$ temperature $\temperature{ik}$ to interface $\interface{ij}$ in \cref{eq:heat_flux_stationary_model_1} indicating that the \textit{stationary} model gives closure law functions which propagate boundary related data \textbf{instantly} between interfaces of the domain.

Let us mention that, in the literature, one can find other nonlinear closure laws like $\heatFlux{ij}(\mathbf{u}_i, \mathbf{b}_{i\star}) \propto \length[\alpha]{i} (\temperature{i}-\temperature{ij})^{\beta}$ with $\alpha,\beta \in \mathbb{R}$ for convective heat transfer or $\heatFlux{ij}(\mathbf{u}_i, \mathbf{b}_{i\star}) \propto {\temperature{ij}}^4$ for radiative heat transfer. In this paper, we will only consider \cref{eq:heat_flux_stationary_model_1,eq:heat_flux_stationary_model_2}.

\section{Introduction}
\label{sec:introduction}
\input{introduction/main.tex}
\section{Governing equations}
\label{sec:governing_equations}
\input{equations/main.tex}
\section{Coupling formulation and schemes}
\label{sec:coupling_formulation}
\input{coupling_scheme/main.tex}
\section{Numerical analysis and examples}
\label{sec:numerical_analysis}
\input{numerical_analysis/main.tex}
\section{Conclusion and perspectives}
\label{sec:conclusion}
\input{conclusion/main.tex}
\section{Acknowledgments}
This work has been carried out within the framework of the \texttt{PROCOR}{} platform development funded by CEA, EDF and AREVA.
\bibliographystyle{phdbib}

\subsection{Linear stability of the coupling of lumped parameter solvers}
\label{sec:discontinuity_numerical_analysis}
First, we consider the Dirichlet-Neumann boundary conditions at the fixed interface $\interface{12}$ given by \cref{eq:thermal_equilibrium_condition_1,eq:thermal_equilibrium_condition_2,eq:thermal_equilibrium_condition_3} $\mFlowRate{12} = -\mFlowRate{21} = 0, \temperature{12} = \temperature{21}, \heatFlux{12} = -\heatFlux{21}$. The two domain masses are constant and only their energy conservation equations are calculated. In the following, we analyze the linear stability of the coupling of LP{} solvers. They are based on a time discretization of the energy conservation \cref{eq:0D_temperature_equation} and closure laws \cref{eq:heat_flux_stationary_model_1,eq:heat_flux_stationary_model_2} to solve the heat conduction problem. A similar stability analysis for finer 1D solvers can be found in \cite{giles_stability_1997}.
\subsubsection{Stability of a \textit{toy} explicit coupling scheme}
\label{sec:toy_explicit_coupling}
The aim here is to build a prototype of a `toy' ECS{} with no subcycling and only one level of time discretization. Such coupling scheme should allow us to solve the coupled problem while solvers of domains $\domain{1}$ and $\domain{2}$ are only called once during a time step. To do so, we use an implicit Euler scheme for the first domain $\domain{1}$ and an explicit Euler scheme for the second domain $\domain{2}$ . Denoting by $\Delta t$ the time step, the discretized equations for the first domain $\domain{1}$ are given by
\begin{equation}
\label{eq:explicit_coupling_scheme_1}
\left\{
	\begin{aligned}
		\heatFlux[n]{12} &= \thermalCond{1} \frac{6 \temperature[n]{1} - 2 \temperature{b_1} - 4 \temperature[n]{12}}{\length{1}}\quad\text{in}\;\domain{1}\\
		\massDensity{1} \heatCp{1} \length{1} \frac{\temperature[n]{1} - \temperature[n-1]{1}}{\Delta t} &= -\heatFlux[n]{12}\quad\text{in}\;\domain{1}
	\end{aligned}
\right.
\end{equation}
with continuity of temperature at the interface, i.e.
\begin{equation}
\label{eq:explicit_coupling_scheme_12}
	\temperature[n]{12} = \temperature[n]{21}\quad\text{on}\;\interface{12}
\end{equation}
and for the second domain $\domain{2}$
\begin{equation}
\label{eq:explicit_coupling_scheme_2}
\left\{
	\begin{aligned}
		\temperature[n+1]{21} &= -\frac{1}{4}\frac{\length{2}}{\thermalCond{2}} \heatFlux[n]{21} + \frac{3}{2} \temperature[n+1]{2} - \frac{1}{2} \temperature{b_2}\quad\text{in}\;\domain{2}\\
		\massDensity{2} \heatCp{2} \length{2} \frac{\temperature[n+1]{2} - \temperature[n]{2}}{\Delta t} &= -\heatFlux[n]{21}\quad\text{in}\;\domain{2}\\
	\end{aligned}
\right.
\end{equation}
with continuity of heat flux at the interface, i.e.
\begin{equation}
\label{eq:explicit_coupling_scheme_21}
		\heatFlux[n]{21} = -\heatFlux[n]{12}\quad\text{on}\;\interface{12}.
\end{equation}
The ECS{} is sequenced as follows : given an interface temperature $\temperature[n]{21}$ at time $t^n$, the first domain is advanced from time $t^{n-1}$ to $t^n$ and a heat flux $\heatFlux[n]{12}$ is computed at time $t^n$ which is then imposed to the second domain. It is then advanced from time $t^n$ to $t^{n+1}$ and it computes a new temperature $\temperature[n+1]{21}$ at time $t^{n+1}$, and so on. Thus, the coupling scheme only requires one call to each solver per time step so that the scheme can be called explicit.

We now analyze the linear stability of the coupling scheme. Combining \cref{eq:explicit_coupling_scheme_1,eq:explicit_coupling_scheme_12,eq:explicit_coupling_scheme_2,eq:explicit_coupling_scheme_21}, we get
\begin{equation}
\label{eq:linear_difference_equation}
	\left(1+6\frac{\Delta t}{\tau_1}\right)\heatFlux[n+1]{12} - \left(1+\left(1-6\frac{\Delta t}{\tau_2}\right)\hbar\right)\heatFlux[n]{12} + \hbar\,\heatFlux[n-1]{12} = 0
\end{equation}
with $\hbar=\frac{\thermalCond{1} / \length{1}}{\thermalCond{2} / \length{2}}$, the characteristic times of conduction $\tau_1 = \massDensity{1} \heatCp{1} {\length{1}}^2 / \thermalCond{1}$ in domain $\domain{1}$ and $\tau_2 = \massDensity{2} \heatCp{2} {\length{2}}^2 / \thermalCond{2}$ in domain $\domain{2}$. To ensure stability of the explicit scheme in domain $\domain{2}$, $\Delta t$ has to be small in comparison to $\tau_2$, i.e.\ $\Delta t / \tau_2 \ll 1$. The simplified characteristic polynomial $\chi$ associated to the simplified second order linear difference \cref{eq:linear_difference_equation} has two complex roots $x_{\{1,2\}}^{\star}$ which have to be of modulus strictly smaller than one in order to guaranty a stable coupling scheme.

We get the following results :
\begin{itemize}
\item When $\Delta t / \tau_1 \ll 1$ and $\hbar<1$, the two roots are real and a Taylor expansion gives
\begin{align*}
	x_1^{\star} &= 1-\frac{6}{1-\hbar}\frac{\Delta t}{\tau_1} + \mathcal{O}(\left(\frac{\Delta t}{\tau_1}\right)^2)\\
	x_2^{\star} &= \hbar+\frac{6\hbar^2}{1-\hbar}\frac{\Delta t}{\tau_1} + \mathcal{O}(\left(\frac{\Delta t}{\tau_1}\right)^2)
\end{align*}
When $\Delta t / \tau_1 \ll 1$ and $\hbar>1$, the two roots are complex of modulus \[\abs{x_{\{1,2\}}^{\star}}^2 = \frac{\hbar}{1+6\frac{\Delta t}{\tau_1}}\]
\item For any values of $\Delta t / \tau_1$, when $\hbar \ll 1$, roots of the characteristic polynomial $\chi$ are $x_1^{\star} \approx \frac{1}{1+6\frac{\Delta t}{\tau_1}}$ and $x_2^{\star} \approx \frac{\hbar}{1+6\frac{\Delta t}{\tau_1}}$ and when $\hbar \gg 1$,  $x_1^{\star} \approx 1+\frac{6\frac{\Delta t}{\tau_1}}{\hbar}$ and $x_2^{\star} \approx \frac{\hbar}{1+6\frac{\Delta t}{\tau_1}}$.
\end{itemize}

Thus, in both cases when $\hbar>1$ the spurious solution for $\heatFlux{12}$ of \cref{eq:linear_difference_equation} associated to root $x_2^{\star}$ increases with a growth rate of $\frac{\hbar}{1+6\,\Delta t / \tau_1}$. In most cases, the spurious solution will grow in time leading to an unstable scheme. On the contrary, cases when $\hbar<1$ give a stable coupling scheme.

Thus, the remaining question is the value $\hbar^{\textrm{crit}}$ of $\hbar$ above which the instabilities appear, thus defining the stability limit region of the ECS{}. Further calculations show that $\hbar^{\textrm{crit}} = \abs{\frac{1+\ 6\ \Delta t/\tau_1}{1-\ 6\ \Delta t / \tau_2}}$ allowing to define a \textit{pseudo} CFL condition for the linear stability of the ECS{} :
\begin{equation}
\label{eq:CFL_explicit}
r_{12}(\Delta t) \overset{\textit{\tiny def.}}{=} \abs{\frac{1-\ 6\ \Delta t / \tau_2}{1+\ 6\ \Delta t/\tau_1}} \hbar < 1
\end{equation}

Finally, the Dirichlet-Neumann explicit coupling should be done in such a way that the domain with the Dirichlet boundary condition, i.e.\ imposed boundary temperature, has the lower thermal conductivity or the higher characteristic length and the domain with the Neumann boundary condition, i.e.\ imposed boundary heat flux, has the higher thermal conductivity and the lower characteristic length. Thus, the stability limit region $\hbar^{\textrm{crit}}$ can be extended for larger values of $\Delta t / \tau_1$ which can be explained by the diffusive properties of the ICS{} used in domain $\domain{1}$.
\paragraph{Numerical evidence of stability.} This analysis is illustrated by \cref{fig:explicit_thermal_interface}. At $t=0^-$, both domains are at an equilibrium temperature of $\temperature{1} = \temperature{2} = \temperature{21} = \temperature{b_1} = \temperature{b_2} = 2000$ K. Boundary discontinuities are imposed at $t=0^+$, $\temperature{b_1} = 3000$ K and $\temperature{b_2} = 400$ K and at $t=3\,\tau_1$, $\temperature{b_1} = \temperature{b_2} = 2000$ K. In addition, all of the computations use the value $\Delta t / \tau_2 = 1/100$ and $\Delta t / \tau_1 = 1/10$ corresponding to a stability limit region $\hbar^{\textrm{crit}} = 1.6$. \Cref{fig:explicit_thermal_interface_11} shows that when $\hbar = \hbar^{\textrm{crit}}$ the oscillations created by the transient initiated by the boundary discontinuities are not damped but are stable while in contrary in \cref{fig:explicit_thermal_interface_21} with $\hbar < \hbar^{\textrm{crit}}$ the oscillations are well damped and in \cref{fig:explicit_thermal_interface_22} with $\hbar > \hbar^{\textrm{crit}}$ the oscillations are clearly unstable. However, those three figures show that the discontinuities at boundaries are instantly propagated to the interface temperature and to the coupled domain in only one coupling time step. For comparison purpose, solution given by the coupling of the 1D heat equation solvers is given in \cref{fig:explicit_thermal_interface_12}. For the computation, the same physical values as for the coupling of LP{} solvers are used. The values meet the requirements expressed in \cite{giles_stability_1997} to ensure proper stability of such solvers. For this coupling, the discontinuities at boundaries start a much slower and smaller transient ($\temperature{21}(t=1000\,\text{s}) \approx 2150$ K) than for the LP{} coupling ($\temperature{21}(t=100\,\text{s}) \approx 2700$ K). The fact that LP{} models propagate in one coupling time step all their discontinuities to the other models tend to create fast and important transient. Thus, the use of ECS{} might not be adapted for this case.

\begin{figure}
    \centering
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_1_a_0_1_b_0_0_1.eps}
		\caption{LP{} solvers, $\hbar=1.6$}
		\label{fig:explicit_thermal_interface_11}
    \end{subfigure}%
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_1D_r_0_2_a_0_0_1_b_0_0_1.eps}
		\caption{multidimensional solvers, $\hbar=1.6$}
		\label{fig:explicit_thermal_interface_12}
    \end{subfigure}\\
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_1_0_2_a_0_1_b_0_0_1.eps}
		\caption{$\hbar=1.63$}
		\label{fig:explicit_thermal_interface_21}
    \end{subfigure}%
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_0_6_2_5_a_0_1_b_0_0_1.eps}
		\caption{$\hbar=1.0$}
		\label{fig:explicit_thermal_interface_22}
    \end{subfigure}
	\caption{Propagation of discontinuities at boundaries in the explicit coupling of lumper parameter and finer (upper right) solvers for heat conduction for different values of $\hbar$ with $\Delta t / \tau_2 = 1/100$ and $\Delta t / \tau_1 = 1/10$.}
	\label{fig:explicit_thermal_interface}
\end{figure}

\subsubsection{Stability of a \textit{toy} implicit coupling scheme}
The iterative algorithm is described below. Given an initial boundary temperature $\temperature[n+1,0]{21}$ for domain $\domain{1}$, it iterates over $k \geq 0$ with the following steps~:
\begin{enumerate}
\item Use an \textbf{implicit} solver in $\domain{1}$ to find the interface heat flux $\heatFlux[n+1,k+1]{12}$
\begin{equation}
\label{eq:implicit_coupling_scheme_1}
\left\{
	\begin{aligned}
		\heatFlux[n+1,k+1]{12} &= \thermalCond{1} \frac{6 \temperature[n+1]{1} - 2 \temperature{b_1} - 4 \temperature[n+1]{12}}{\length{1}}\quad\text{in}\;\domain{1}\\
		\massDensity{1} \heatCp{1} \length{1} \frac{\temperature[n+1]{1} - \temperature[n]{1}}{\Delta t} &= -\heatFlux[n+1,k+1]{12}\quad\text{in}\;\domain{1}
	\end{aligned}
\right.
\end{equation}
with continuity of temperature at the interface, i.e.
\begin{equation}
\label{eq:implicit_coupling_scheme_12}
	\temperature[n+1]{12} = \temperature[n+1,k]{21}\quad\text{on}\;\interface{12}.
\end{equation}
\item Use an \textbf{implicit} solver in $\domain{2}$ to find the interface temperature $\tilde{T}^{n+1,k+1}_{21}$%
\begin{equation}
\label{eq:implicit_coupling_scheme_2}
\left\{
	\begin{aligned}
		 \tilde{T}^{n+1,k+1}_{21} &= -\frac{1}{4}\frac{\length{2}}{\thermalCond{2}} \heatFlux[n+1]{21} + \frac{3}{2} \temperature[n+1]{2} - \frac{1}{2} \temperature{b_2}\quad\text{in}\;\domain{2}\\
		 \massDensity{2} \heatCp{2} \length{2} \frac{\temperature[n+1]{2} - \temperature[n]{2}}{\Delta t} &= -\heatFlux[n+1]{21}\quad\text{in}\;\domain{2}\\
	\end{aligned}
\right.
\end{equation}
with continuity of heat flux at the interface, i.e.
\begin{equation}
\label{eq:implicit_coupling_scheme_21}
		\heatFlux[n+1]{21} = -\heatFlux[n+1,k+1]{12}\quad\text{on}\;\interface{12}.
\end{equation}
\item Measure convergence level with
\begin{equation}
	\frac{\abs{\tilde{T}^{n+1,k+1}_{21} - \temperature[n+1,k]{21}}}{\temperature[n+1,k]{21}} \leq \epsilon_{\textrm{rel}}.
\end{equation}
If the previous predicate is not satisfied, relax the interface temperature by
\begin{equation}
	\temperature[n+1,k+1]{21} = \omega \tilde{T}^{n+1,k+1}_{21} + (1-\omega)\temperature[n+1,k]{21}
\end{equation}
and loop again. Otherwise, the final interface heat flux and temperature are given by $\temperature[n+1]{21} = \temperature[n+1,k]{21}$ and $\heatFlux[n+1]{12} = \heatFlux[n+1,k]{12}$. Consequently, internal variables $\temperature[n+1]{1}$ and $\temperature[n+1]{2}$ can be fully implicitly calculated : the iterative algorithm effectively allows to use implicit solvers for the two domains while allowing to decouple the two domains.
\end{enumerate}

Combining \cref{eq:implicit_coupling_scheme_1,eq:implicit_coupling_scheme_12,eq:implicit_coupling_scheme_2,eq:implicit_coupling_scheme_21}, simple calculations lead to the following equation giving the interface temperature iterates $\temperature[n+1,k]{21}$
\begin{equation}
\temperature[n+1,k+1]{21} = \left[1-\left(1+\frac{1-6\frac{\Delta t}{\tau_2}}{1+6\frac{\Delta t}{\tau_1}}\hbar\right)\omega\right]\temperature[n+1,k]{21} + g(\temperature[n]{1}, \temperature[n]{2}, \temperature{b_1}, \temperature{b_2}).
\end{equation}
We note $\temperature[n+1]{21}$ the solution of the fixed point equation associated with the previous equation and $e_k =\abs{\temperature[n+1,k]{21} - \temperature[n+1]{21}}$ the error at iteration $k$. Then iteration errors are given by equation
\begin{equation}
\label{eq:error_toy_implicit_coupling}
e_{k+1} = \abs{1-\left(1+\frac{1-6\frac{\Delta t}{\tau_2}}{1+6\frac{\Delta t}{\tau_1}}\hbar\right)\omega} e_k.
\end{equation}
Thus, the iterative algorithm converges toward the interface solutions, i.e.\ $\temperature[n+1,\infty]{21} = \temperature[n+1]{21}$ and $\heatFlux[n+1, \infty]{12} = \heatFlux[n+1]{12}$, if and only if
\begin{equation}
\label{eq:omega_condition_convergence}
0 < \omega < \frac{2}{1+\frac{\abs{1-6\frac{\Delta t}{\tau_2}}}{\abs{1+6\frac{\Delta t}{\tau_1}}}\hbar} = \frac{2}{1+r_{12}(\Delta t)}.
\end{equation}
Note that when $\frac{\Delta t}{\tau_1} \ll 1$ and  $\frac{\Delta t}{\tau_2} \ll 1$, condition \cref{eq:omega_condition_convergence} reads
\begin{equation}
	0 < \omega < \frac{2}{1+\hbar}.
\end{equation}
In cases where $\hbar>1$ for which the ECS{} diverges, the iterative algorithm needs $\omega < 1$, i.e. under-relaxation to converge, which might lead to slow convergence.

\paragraph{Numerical evidence of stability.} This analysis is confirmed by \cref{fig:implicit_thermal_interface}. As in \cref{sec:toy_explicit_coupling}, both domains are at an equilibrium temperature of $2000$ K and the same discontinuities at boundaries are occurring at the same time $t=0^{+}$ and $t=3\,\tau_1$. For the computations, we still use $\Delta t / \tau_2 = 1/100$, $\Delta t / \tau_1 = 1/10$. We use $\hbar=1.6$ corresponding to the stability limit region for the ECS{} ($\hbar = \hbar^{\textrm{crit}} = 1.6$) for which this scheme is showing constant oscillations. For this value, the corresponding maximal value of $\omega$ to ensure convergence of the iterative process is $\omega \approx 1.03$. The relative tolerance used is $\epsilon_{\textrm{rel}}=10^{-4}$. Several computations confirm the predicted behavior of the ICS{} : as long as $\omega$ stays under the calculated maximum value $1.03$, the ICS{} converges toward the same smooth solution given by \cref{fig:implicit_thermal_interface_11}. Even for values $\hbar > \hbar^{\textrm{crit}}$, the ICS{} handles smoothly the fast dynamics and important transient initiated by the instant propagation of the discontinuities at boundaries in the two domains by the LP{} models. For cases where the ECS{} is not unconditionally unstable, i.e. $r_{12}(\Delta t) < 1$, one way to achieve stability is to reduce the coupling time step : \cref{fig:implicit_thermal_interface_12} shows an explicit coupling solution with $\hbar = 0.625 < \hbar^{\textrm{crit}} \approx 1$ and coupling time step $\Delta t$ reduced by a factor of $10$. However, while effectively reducing the oscillations, using an ECS{} with such a small coupling time step leads to higher computational time. The ICS{} allows to keep a sufficiently large coupling time step and a limited number of iterations as shown in \cref{fig:implicit_thermal_interface_21,fig:implicit_thermal_interface_22}. In the last figure, it is important to understand that only the solution given by the ECS{} with a small coupling time step is usable but is costly to obtain while the solution given with the ICS{} is good and costs $5$ times less. For this problem, the optimal relaxation parameter leading to the higher convergence rate can be analytically calculated from \cref{eq:error_toy_implicit_coupling} and is given by $\omega = 1/(1+r_{12}(\Delta t)) \approx 0.52$. However, most of the time the optimal relaxation parameter is highly dependent to the problem.

\begin{figure}
    \centering
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_implicit_r_1_a_0_1_b_0_0_1.eps}
		\caption{ICS{}, $\hbar=1.6$}
		\label{fig:implicit_thermal_interface_11}
    \end{subfigure}%
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_0_6_2_5_a_0_0_1_b_0_0_1.eps}
		\caption{ECS{}, $\hbar=0.625$, $\Delta t \leftarrow \Delta t / 10$}
		\label{fig:implicit_thermal_interface_12}
    \end{subfigure}\\
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_1_a_0_1_b_0_0_1_iterations.eps}
		\caption{ICS{} number of iterations, $\hbar=1.6$}
		\label{fig:implicit_thermal_interface_21}
    \end{subfigure}%
	\begin{subfigure}[b]{.5\textwidth}
        \includegraphics[width=\textwidth]{thermal_interface_r_1_a_0_1_b_0_0_1_simulation_times.eps}
		\caption{Schemes time of simulation, $\hbar=1.6$}
		\label{fig:implicit_thermal_interface_22}
    \end{subfigure}
	\caption{Propagation of discontinuities at boundaries of lumper parameter heat conduction coupled solvers solved with an ICS{} (upper left) and an ECS{} (upper right) with much smaller coupling time step.}
	\label{fig:implicit_thermal_interface}
\end{figure}

\subsection{Synchronization on internal events}
\label{sec:synchronisation_numerical_analysis}
For instance, let us describe the state of domain $\domain{2}$ into a graph of three states as depicted in \cref{fig:states_toy_model} : \textit{Heating}, \textit{Melting} and \textit{Empty} states. A state transition occurs when a certain algebraic function is activated, e.g. $\mFlowRate{21} < 0$. As explained \cref{sec:coupling_scheme}, each state has its own solver, its own set of equations and its own interface and boundary conditions. State transitions can lead to discontinuities of state or interface variables.
\begin{figure}
	\centering
	\begin{tikzpicture}[node distance=4.5cm]
		\tikzstyle{vertex}=[ellipse, thick, draw=black, text width=1.5cm, align=center]
		\node[vertex, fill=red!50](h) at (0,0) {$\mathcal{M}^{\text{Heating}}$};
 		\node[vertex, fill=red!70, right of=h](m) {$\mathcal{M}^{\text{Melting}}$};
 		\node[vertex, fill=gray!10, right of=m](e) {$\mathcal{M}^{\text{Empty}}$};
		\path (h) edge[->, thick, bend right=10, dashed] node[midway, below] {$\temperature{21} \geq \temperature[\star]{}$} (m);
		\path (m) edge[->, thick, bend right=10, dashed] node[midway, above] {$\mFlowRate{21} < 0$} (h)
				  edge[->, thick, dashed] node[midway, below] {$\mass{2} \leq \mass[\epsilon]{} \approx 0$} (e);
	\end{tikzpicture}
	\caption{Internal states of domain $\domain{2}$ and their transition functions.}
	\label{fig:states_toy_model}
\end{figure}
For this case, equations of solver $\mathcal{M}^{\textrm{Heating}}$ are the ones described previously in \cref{sec:discontinuity_numerical_analysis} while solver $\mathcal{M}^{\textrm{Melting}}$ continuous equations are given by
\begin{equation}
\left\{
	\begin{aligned}
		\heatFlux{21} &= \thermalCond{1} \frac{6 \temperature{2} - 4 \temperature{21} - 2 \temperature{b_2}}{\length{2}}\quad\text{in}\;\domain{2}\\
		\D{\mass{2}}{t} &= \massDensity{2}\D{\length{2}}{t} = -\mFlowRate{21}	\quad\text{in}\;\domain{2}\\
		\massDensity{2} \heatCp{2} \length{2} \D{\temperature{2}}{t} &= -\heatFlux{21}\quad\text{in}\;\domain{2},\\
	\end{aligned}
\right.
\end{equation}
leading to a mobile fusion solidification front at interface $\interface{12}$ given by \cref{eq:fusion_front_equilibrium_condition_1,eq:fusion_front_equilibrium_condition_2,eq:fusion_front_equilibrium_condition_3} $\heatFlux{12} = - \heatFlux{21} + \dEnthalpy[\textrm{fus.}]{} \mFlowRate{21}$, $\mFlowRate{12} = -\mFlowRate{21}$ and $\temperature{12} = \temperature{21} = \temperature[\textrm{fus.}]{}$. Again, another set of internal and boundary equations are used for solver $\mathcal{M}^{\textrm{Empty}}$ leading to a fixed interface given by \cref{eq:thermal_equilibrium_condition_1,eq:thermal_equilibrium_condition_2,eq:thermal_equilibrium_condition_3} $\mFlowRate{12} = -\mFlowRate{21} = 0, \temperature{12} = \temperature{21}, \heatFlux{12} = -\heatFlux{21}$.

Again, at $t=0^-$, both domains are at an equilibrium temperature of $\temperature{1} = \temperature{2} = \temperature{21} = \temperature{b_1} = \temperature{b_2} = 2000$ K. In particular, domain $\domain{2}$ is in \textit{Heating} state. At $t=0^+$, $\temperature{b_1} = 3000$ K and $\temperature{b_2} = 3000$ K to force the \textit{Melting} state transition when $\temperature{21}$ reaches $\temperature[\star]{} = 2100$ K at time $t_{H \rightarrow M}^{\star}$. Once this state has been reached, the second domain is in \textit{Melting} state until only a residual mass $\mass{2} = \mass[\epsilon]{} = 150$ kg is left (corresponding to height $\length{2} = 1.5$ cm) and the \textit{Empty} state is reached at time $t_{M \rightarrow E}^{\star}$. The computation then stops when domain $\domain{1}$ reaches a near stationary state. In addition, all of the computations use the value $\tau_2 = 10^4$ s, $\tau_1 = 8\times10^3$ s and different values for the time step $\Delta t$. The chosen values ensure that the ECS{} is stable, i.e. $r_{12}(\Delta t) < 1.0$. \Cref{tab:explicit_implicit_synchro} shows the times of the internal events of domain $\domain{2}$ computed by the ECS{} and the ICS{} for several coupling time steps. Usually, we ask the coupling time step to be bounded by \[\tau / 10 = 800\ \text{s} \geq \Delta t \geq \tau/100 = 80\ \text{s}\quad\text{with}\quad\tau=\min{}(\tau_1,\tau_2)\] for computational efficiency (in red and blue in the table). In comparison to the reference solution (marked with $\dagger$) given by the ECS{} for a small time step $\Delta t= 1$ s for which the scheme has converged, the ICS{} is able to predict the state transition with high accuracy and adapt its time step to synchronize both domains when the event occurs. However, the synchronization mechanism of the ECS{} only allows the domains to be informed of an event at the end of the time step. During this window of time, both domains are in a non physical state, i.e.\ the domain $\domain{1}$ is not aware of the disappearance of the domain $\domain{2}$. This leads to numerical errors of order $\mathcal{O}(\Delta t)$ and the numerical creation of either mass and/or energy leading to different stationary states for the ECS{}. For a time step $\Delta t = 100$ s respecting the constraint, the explicit solution leads to a relative error of $18 \%$ in term of mass. Reducing the coupling time step to $10$ s leads to a relative error of $0.3 \%$ but, as seen previously, increases computational time. Besides, the ICS{} always converges toward the reference solution and the right stationary state even for large coupling time step.\\
\begin{table}
\centering
{\renewcommand{\arraystretch}{1.2}
\begin{tabular}{l|c|c|c|c|c|c}
& $\Delta t$ (s) & $\textcolor{blue}{100}$ & $\textcolor{purple}{50}$ & $\textcolor{purple}{25}$ & $\textcolor{purple}{10}$ & $\textcolor{purple}{1}$\\
\hline
& $r_{12}(\Delta t)$ & \multicolumn{5}{c}{< 1.0}\\
\hline
\hline
\multirow{2}{*}{ECS{}} & $\mass[\infty]{1}$ (kg) & $\textcolor{blue}{1492}$ & $\textcolor{purple}{1413}$& $\textcolor{purple}{1347}$ & $\textcolor{purple}{1260}$ & $\textcolor{purple}{1256}\ \dagger$\\
\cline{2-7}
& $\temperature[\infty]{1}$ (K) & $\textcolor{blue}{1710}$ & $\textcolor{purple}{1708}$& $\textcolor{purple}{1705}$ & $\textcolor{purple}{1702}$ & $\textcolor{purple}{1700}\ \dagger$\\
\hline
\multirow{2}{*}{ICS{}} & $\mass[\infty]{1}$ (kg) & $\textcolor{blue}{1255}$ &\multicolumn{4}{c}{$\textcolor{purple}{1255}$}\\
\cline{2-7}
& $\temperature[\infty]{1}$ (K) & $\textcolor{blue}{1700}$ & \multicolumn{4}{c}{$\textcolor{purple}{1700}$}\\
\hline
\hline
\multirow{2}{*}{ECS{}} & $t_{H \rightarrow M}^{\star}$ (s) & $\textcolor{blue}{1700}$ & $\textcolor{purple}{1600}$& $\textcolor{purple}{1600}$ & $\textcolor{purple}{1590}$ & $\textcolor{purple}{1584}\ \dagger$\\
\cline{2-7}
& $t_{M \rightarrow E}^{\star}$ (s) & $\textcolor{blue}{3200}$ & $\textcolor{purple}{3200}$& $\textcolor{purple}{3150}$ & $\textcolor{purple}{3120}$ & $\textcolor{purple}{3117}\ \dagger$\\
\hline
\multirow{2}{*}{ICS{}} & $t_{H \rightarrow M}^{\star}$ (s) & $\textcolor{blue}{1583}$ & \multicolumn{4}{c}{$\textcolor{purple}{1583}$}\\
\cline{2-7}
& $t_{M \rightarrow E}^{\star}$ (s) & $\textcolor{blue}{3120}$ & $\textcolor{purple}{3117}$& $\textcolor{purple}{3117}$ & $\textcolor{purple}{3116}$ & $\textcolor{purple}{3116}$\\
\hline
\end{tabular}}
\caption{Times $t_{H \rightarrow M}^{\star}$ and $t_{M \rightarrow E}^{\star}$ of domain $\domain{2}$ internal events and stationary state $\mass[\infty]{1}$ and $\temperature[\infty]{1}$ reached by domain $\domain{1}$. Reference solution marked with $\dagger$, solutions respecting the coupling time step constraint in \textcolor{blue}{blue}, otherwise in \textcolor{purple}{red}.}
\label{tab:explicit_implicit_synchro}
\end{table}

This very simple example of coupling of LP{} models highlights the main drawback of such lightweight modeling : deleting spatial dependency in the equations force each model to instantly propagate all its data, e.g.\ its discontinuities or its state transitions, thus creating high and important transient or bringing the system in non physical states. Fixes are used to bring back the system in a coherent state which should be avoided at all cost. These phenomenons create numerical errors in the coupling which can be measured in term of numerical energy created at the interface $\interface{12}$. We follow the same methodology as in \cite{piperno_partitioned_2001}. The local variation of energy $\Delta E^{n \rightarrow n+1}$ at interface $\interface{12}$ between times $t^n$ and $t^{n+1}$ is the sum of the energy $\Delta E_1^{n \rightarrow n+1}$ send by domain $\domain{1}$ through the interface and the energy $\Delta E_2^{n \rightarrow n+1}$ send by domain $\domain{2}$ through the interface. For the continuous case in which both domains are strongly coupled, the local variation of energy is null. It is defined by equation
\begin{align*}
\Delta E^{n \rightarrow n+1} &= \Delta E_1^{n \rightarrow n+1}\ -\ (-\Delta E_2^{n \rightarrow n+1})\\
&=\int_{t^n}^{t^{n+1}} \heatFlux{12} - \left(\int_{t^n}^{t^{n+1}} \heatFlux{21} - \dEnthalpy[\textrm{fus.}]{} \mFlowRate{21}\right)\\
&= 0.
\end{align*}
The global energy $\Delta E^{0 \rightarrow n}$ through interface $\interface{12}$ is defined by
\begin{equation}
\label{eq:global_energy}
\Delta E^{0 \rightarrow n} = \sum_{k=0}^{n} \Delta E^{k \rightarrow k+1} = 0
\end{equation}
and is also null in the continuous case. For each energy variation $\Delta E$, we define the relative energy
\begin{equation*}
\epsilon(\Delta E) = \frac{\Delta E}{E^{\star}} = \frac{\Delta E}{\mass{2}(t^0)\heatCp{2}(\temperature[\textrm{fus.}]{} - \temperature{2}(t^0)) + \dEnthalpy[\textrm{fus.}]{}\mass{2}(t^0)}
\end{equation*}
by the ratio of the energy $\Delta E$ to the energy $E^{\star}$ needed to heat the initial domain $\domain{2}$ to its fusion temperature ($\mass{2}(t^0)\heatCp{2}(\temperature[\textrm{fus.}]{} - \temperature{2}(t^0))$) and to melt it ($\dEnthalpy[\textrm{fus.}]{}\mass{2}(t^0)$). With the physical values used, we have $E^{\star} \approx 10^8$ J.

When discretized for the ECS{}, the local variation of energy becomes
\begin{equation*}
\Delta E^{n \rightarrow n+1} \approx \heatFlux{12}(t^{n+1})\Delta t\ -\ \left(\heatFlux{21}(t^n)-\dEnthalpy[\textrm{fus.}]{} \mFlowRate{21}(t^n)\right) \Delta t
\end{equation*}
which emphasizes the ``time-lag'' between the two domains. \Cref{fig:explicit_interface_energy_11} shows that the ECS{} can create locally $2.5 \%$ of $E^{\star}$, i.e. $2.5\,10^6$~J. \Cref{fig:explicit_interface_energy_12} shows first that the global energy is not null and that terms of the sum \cref{eq:global_energy} do not sum to zero. Besides, the explicit scheme can create globally more than $6 \%$ of $E^{\star}$, i.e. $6\,10^6$~J which can have disastrous effects on the whole accuracy.
\begin{figure}
    \centering
	\begin{subfigure}[b]{0.5\textwidth}
    	\includegraphics[width=1.0\textwidth]{explicit_interface_local_energy.eps}
		\caption{}
		\label{fig:explicit_interface_energy_11}
    \end{subfigure}%
	\begin{subfigure}[b]{0.5\textwidth}
    	\includegraphics[width=1.0\textwidth]{explicit_interface_global_energy.eps}
		\caption{}
		\label{fig:explicit_interface_energy_12}
    \end{subfigure}
	\caption{Relative local energy variation (\textit{left}) and global energy variation (\textit{right}) at interface for the ECS{}}
	\label{fig:explicit_interface_energy}
\end{figure}

For the ICS{}, the local variation of energy becomes
\begin{equation*}
\Delta E^{n \rightarrow n+1} \approx \heatFlux{12}(t^{n+1,\infty})\Delta t\ -\ \left(\heatFlux{21}(t^{n+1,\infty})-\dEnthalpy[\textrm{fus.}]{} \mFlowRate{21}(t^{n+1,\infty})\right) \Delta t
\end{equation*}
in which variables with superscript $^\infty$ represent the value reached when the convergence criterion of the iterative scheme is satisfied. At convergence, we can bound the residual at interface in such a way that
\begin{equation*}
\abs{\heatFlux{12}(t^{n+1,\infty})\ -\ \left(\heatFlux{21}(t^{n+1,\infty})-\dEnthalpy[\textrm{fus.}]{} \mFlowRate{21}(t^{n+1,\infty})\right)}\Delta t \leq \epsilon_{\textrm{rel}}\abs{\heatFlux{12}(t^{n+1,\infty})}\Delta t.
\end{equation*}
This allows us to bound the relative local energy at interface of the ICS{} with
\begin{equation*}
\abs{\epsilon(\Delta E^{n \rightarrow n+1})} = \abs{\frac{\Delta E^{n \rightarrow n+1}}{E^{\star}}} \leq \epsilon_{\textrm{rel}}\ \frac{\abs{\heatFlux{12}(t^{n+1,\infty})}\Delta t}{E^{\star}}
\end{equation*}
and the relative global energy with
\begin{equation*}
\abs{\epsilon(\Delta E^{0 \rightarrow n})} \leq \sum_{k=0}^{n} \abs{\frac{\Delta E^{k \rightarrow k+1}}{E^{\star}}} \leq n\ \epsilon_{\textrm{rel}}\ \frac{\max\limits_{k}\abs{\heatFlux{12}(t^{k+1,\infty})}\Delta t}{E^{\star}}.
\end{equation*}
This shows that the imbalance of energy at interface of the ICS{} can be controlled with the relative tolerance $\epsilon_{\textrm{rel}}$. This is confirmed by numerical experiments presented \cref{fig:implicit_interface_energy_11} and \cref{fig:implicit_interface_energy_12} showing respectively the relative local and global energy at interface with a relative tolerance set to $10^{-4}$. The maximum heat flux $\max_{k}\abs{\heatFlux{12}(t^{k+1,\infty})}$ reached at the interface with unit area is equal to $\approx 1.5\times10^5\ \text{W.m$^{-2}$}{}$ leading to a bound for the relative local energy at interface of $1.5\times10^{-3}$ (expressed in percent in \cref{fig:implicit_interface_energy_11}).
\begin{figure}
    \centering
	\begin{subfigure}[b]{0.5\textwidth}
    	\includegraphics[width=1.0\textwidth]{implicit_interface_local_energy.eps}
		\caption{}
		\label{fig:implicit_interface_energy_11}
    \end{subfigure}%
	\begin{subfigure}[b]{0.5\textwidth}
    	\includegraphics[width=1.0\textwidth]{implicit_interface_global_energy.eps}
		\caption{}
		\label{fig:implicit_interface_energy_12}
    \end{subfigure}
	\caption{Relative local energy variation (\textit{left}) and global energy variation (\textit{right}) at interface for the ICS{}}
	\label{fig:implicit_interface_energy}
\end{figure}

redpajama_set_name: RedPajamaArXiv


Last night I headed down to the Red Poppy Art House for my break between Siegfried and Götterdämmerung during this week's San Francisco Opera performance of Richard Wagner's Der Ring des Nibelung (the ring of the Nibelung). The show for the evening featured flamenco dancer Melissa Cruz in a program entitled Collaborations in Flamenco. The program was given two performances, and I was there for the first of them.
This was definitely not a typical flamenco evening where dancers performed as part of a group of guitarists and singers. Cruz opened the program with a solo performance, an approach to flamenco that I had never previously encountered. While the rhetoric of her facial expressions and upper body language were familiar, the absence of any other sounds led me to concentrate on her footwork and her capacity to summon up elaborate rhythms that could stand securely as music unto itself. Indeed, listening to the sonic language of her feet, so to speak, took me back some 30 years when I would go into Manhattan to listen to Max Roach fill an entire evening with solo drum work.
Once Cruz had firmly established the voice of her own feet, she began to work with different combinations of four musicians who had joined her for the evening. These were percussionist Marco Peris, Bob Sanders on guitar and sometimes percussion, Gregory Masako Jenkins on clarinet and vocal work, and Alex Conde on piano and one percussion gig. Peris was the only one with a conventional drum. The other musicians played a touch-sensitive electronic device in the form of a cube also serving as a stool.
Conde's music was the only offering that came closest to Hispanic idioms. Both Jenkins and Peris worked from different Balkan backgrounds, while Sanders' technique tended to reflect primarily on jazz techniques from the middle of the last century. What was most important, however, was that every combination of musicians providing backup for Cruz had its own characteristic sense of spontaneity; so that each of the selections (which Cruz called vignettes) allowed for its own approaches to improvisation. The program then concluded with the entire group improvising off the tops of their respective heads with Cruz adding her own rhythms and dance style to the fabric of the emerging polyphony.
In that kind of a setting, it was easy to appreciate that Cruz' skills as a musician were just as impressive as the personal stamp she put on her flamenco dancing.


Happy Friday babes! Today I am sharing with you a machine that allows you to do Airbrush makeup at home! It is called Luminess Air. It is the #1 Airbrush makeup system and it is so easy to use! I was honestly nervous at first because I have never used one before but I read through the tutorial guide in the book provided and watched the online application tutorial and I think I did pretty good! They give you a variety of skin tone shades and you can try them out on your hand before applying them. I loved the fact that is was a no mess application and I used so much less makeup! I always feel like I have to use a lot of foundation but with this you can just use 2-4 drops and it will cover your face!
It blended very evenly for me and I feel like it covered all of my pores and blemishes (YAY)! Luminess Air is known to last up to 18 hours and leave your skin feeling fresh and flawless. I don't know about you guys but I always forget about washing my makeup brushes regularly but with this kit you won't ever have dirty brushes touch your face again!
The Airbrush system I am using is here and it is $50 off right now! Plus the black and white is sooo cute!
Thank you to Luminess Air for sponsoring today's post!
This sounds kind of amazing! You look gorgeous, babe.
i had no idea you could do this at home!!
This sounds awesome! I'm going to check out their website. Love the photos! So gorgeous!
No way thats so cool!!!
This looks so cool! I'll have to check it out!
I've been wanting to try one of these on my own! I've had my makeup airbrushed by someone before but haven't done it myself! How cool!
I've always wanted to try airbrush makeup!
This sounds amazing! I've always wanted to try out airbrush makeup!
I've heard so much about this product! I would love to try it!
This sounds amazing! Would love to try it out sometime!
Great shared!! airbrushes are just mind blowing. I would like to share a product which I had used recently and it had worked great on my skin that's dermology anti aging cream. It's a great product and I would recommend it to every one.


Clymer Manuals Suzuki LT230 and LT250 M475 Suzuki LT230S Chain Drive, LFF230G Shaft Drive, LT230G Shaft Drive and LT250S Chain Drive manual. Clymer ATV repair manuals are written specifically for the doityourself enthusiast. Official Suzuki parts list Easy repairs with OEM diagrams Free acccess to parts fiches for Suzuki LT250EF 1985 (F) 85 Suzuki Lt250ef Atv Service Manual 85 Suzuki Lt250ef Atv Service Manual In this site is not the same as a solution directory you buy in a tape gathering or download off the web.
Our on top of 12, 867 manuals and Ebooks is the reason why customers save If you are searching for a ebook 85 suzuki lt250ef manual in pdf format, then you have come on to faithful site. We present full edition of this book in 85 Suzuki Lt250ef Atv Service Manual 85 Suzuki Lt250ef Atv Service Manual In this site is not the thesame as a answer manual you buy in a compilation gathering or download off the web.
Our over 9, 603 manuals and Ebooks is the excuse why Sep 02, 2018 ATV General Discussion General Discussions about Suzuki ATVs. Suzuki ATV Forum Suzuki ATV Discussions ATV General Discussion Needed 1985 Suzuki LT250EF Manual download link Reply Suzuki LT250E LT250EF 1985 Workshop Service Repair 85 suzuki lt250ef manual Download This is the COMPLETE official full factory service repair manual for the Suzuki LT250E LT250EF 1985.
Hundreds of pages allow you to print it out Suzuki atv manuals Clymer Suzuki manuals are written specifically for the doityourself enthusiast. From basic maintenance and troubleshooting to complete overhauls, our Suzuki manuals provide the information you need. 85 Suzuki Lt250ef Atv Service Manual. pdf Free Download Here TOOLS& SUPPLIES CLYMER MANUALS CanAm ATV Forum Get the parts you need for your 1985 SUZUKI LT250EF QUAD RUNNER at Rocky Mountain ATVMC.
We offer a huge selection of parts and accessories to keep your machine running at peak levels. Details about Suzuki Quick Reference Service Manual Data Sheet 85 LT250EFF LT250EF Quadrunner Be the first to write a review.
Suzuki Quick Reference Service Manual Data Sheet 85 LT250EFF LT250EF Quadrunner Complete service repair workshop manual for the: Suzuki Quad Runner 250 LT250E LT250EF LT250 This is the same manual dealerships use to repair your atv quad.