11institutetext:  Siemens Industry Software, Roanne, France 11email:
<EMAIL_ADDRESS>
22institutetext: Institut Camille Jordan, Université de Lyon, UMR5208
CNRS-U.Lyon1, Villeurbanne, France
22email<EMAIL_ADDRESS>

# IFOSMONDI Co-simulation Algorithm with Jacobian-Free Methods in
PETSc††thanks: Supported by organization Siemens Industry Software.

Yohan Eguillon 1122 0000-0002-9386-4646 Bruno Lacabanne 11 0000-0003-1790-3663
Damien Tromeur-Dervout 22 0000-0002-0118-8100

###### Abstract

IFOSMONDI iterative algorithm for implicit co-simulation of coupled physical
systems (introduced by the authors in july 2019 during the Simultech
conference, p.176-186) enables us to solve the non-linear coupling function
while keeping the smoothness of interfaces without introducing a delay.
Moreover, it automatically adapts the size of the steps between data exchanges
among the systems according to the difficulty of the solving of the coupling
constraint. The latter was solved by a fixed-point algorithm in the original
implementation whereas this paper introduces the JFM version (standing for
Jacobian-Free Methods). Most implementations of Newton-like methods require a
jacobian matrix which can be difficult to compute in the co-simulation
context, except in the case where the interfaces are represented by a Zero-
Order-Hold (ZOH). As far as IFOSMONDI coupling algorithm uses Hermite
interpolation for smoothness enhancement (up to Third-Order-Hold), we propose
hereafter a new formulation of the non-linear coupling function including both
the values and the time-derivatives of the coupling variables. This
formulation is well designed for solving the coupling through jacobian-free
Newton type methods. Consequently, successive function evaluations consist in
multiple simulations of the systems on a co-simulation time-step using
rollback. The orchestrator-workers structure of the algorithm enables us to
combine the PETSc framework on the orchestrator side for the non-linear
Newton-type solvers with the parallel integrations of the systems on the
workers side thanks to MPI processes. Different non-linear methods will be
compared to one another and to the original fixed-point implementation on a
newly proposed 2-systems academic test-case (mass-spring-damper type) with
direct feedthrough on both sides.

###### Keywords:

Co-simulation Systems coupling Coupling methods Jacobian-free Newton PETSc
Parallel integration Strong coupling test case

##  1 Introduction

The use of co-simulation is increasing in the industry as it enables to
connect and simulate systems with given interfaces (input and output
variables) without disclosing the expertise inside. Hence, modellers can
provide system architects with virtual systems as black boxes since the
systems are able to interact through their interfaces. Among these
interactions, the minimal requirement are quite simple: a system should at
least be able to read the inputs given by the other systems, to simulate its
physics inside (most of the time thanks to an embedded solver), and to provide
outputs of the simulation to the other systems.

Besides its black box aspect protecting the know-how, co-simulation also
enables physic-based decomposition (one system can represent the hydraulic
part of a modular model, another the mechanical part, a third one the
electrical part, and so on) and/or dynamics-based decomposition (some systems
isolate the stiff state variables so that they do not constraint all the other
states anymore during the simulation). In other words, the co-simulation opens
many doors thanks to the modular aspect of the models handled.

The co-simulation field of research nowadays focuses on the numerical methods
and algorithms that can be used to process simulations of such modular models.
From the simplest implementations (non-iterative Jacobi) to very advanced
algorithms [7, 11, 12, 14, 4], co-simulation methods have been developped in
different fields, showing that the underlying problems to be tackled are not
straightforward. Some arising problems could clearly be identified since the
moment it has become a center of interest for researchers, such as the delay
between the given inputs and the retrieved outputs of a system (corresponding
to the so-called ”co-simulation step” or ”macro-step”), the instabilities that
might occur as a consequence of this delay [15], the discontinuities produced
at each communication [5], the error estimation (and the use of it in order to
adapt the macro step size) [12], the techniques to solve the so-called
”constraint function” corresponding to the interface of the systems [9, 13],
and so on. Many of these problems have been addressed in papers either
proposing an analysis, a method to solve them, or both.

In our previous paper [6], an iterative method that satisfies the interfaces
consistency while avoiding discontinuities at each macro-step was proposed and
compared to well-established methods (non-iterative Jacobi, zero-order hold
iterative co-simulation [9], and non-iterative algorithm enhancing variables
smoothness [5]). This algorithm was based on a fixed-point iterative method.
Its evolution, presented in this paper, is based on iterative methods that
normally require jacobian matrix computation, yet we use their jacobian-free
version. The name of this method is IFOSMONDI-JFM, standing for Iterative and
Flexible Order, SMOoth and Non-Delayed Interfaces, based on Jacobian-Free
Methods. The enhancements it brings to the classical IFOSMONDI method enable
to solve cases that could not be solved by this previous version. The
integration of an easily modulable jacobian-free method to solve the
constraint function will be presented. In particular, the software integration
was made possible thanks to the PETSc framework, a library that provides
modulable numerical algorithms. The interfacing between PETSc and the co-
simulation framework dealing with the systems, interfaces and polynomial
representations will be detailed.

##  2 Formalism and notations

###  2.1 A word on JFM accronym

In the whole paper, the abbreviation JFM will denote jacobian-free versions of
iterative methods that are designed to bring a given function (so-called
callback) to zero and that normally require the computation of the jacobian
matrix of the callback function. In particular, a fixed-point method does not
meet these criteria: it is not a JFM, contrary to matrix-free versions of the
Newton method, the Anderson method [1] or the non-linear GMRES method [10].

###  2.2 General notations

In this paper, we will focus on the explicit systems. In other words, we will
consider that every system in the co-simulation is a dynamical system
corresponding to an ODE (Ordinary Differential Equation). The time-domain of
the ODEs considered will be written $[t^{\text{init}},t^{\text{end}}[$, and
the variable $t$ will denote the time.

Let’s consider $n_{sys}\in\mathbb{N}^{*}$ systems are involved: we will use
the index $k\in[\\![1,n_{sys}]\\!]$ to denote the $k$th system, and
$n_{st,k}$, $n_{in,k}$, and $n_{out,k}$ will respectively denote the number of
state variables, the number of inputs, and the number of outputs of system
$k$.

The time-dependant vectors of states, inputs and outputs of system $k$ will
respectively be written $x_{k}\in
L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{st,k}})$, $u_{k}\in
L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{in,k}})$, and $y_{k}\in
L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{out,k}})$ where $L(A,B)$
denotes the set of functions of domain $A$ and co-domain $B$. Finally, we can
write the ODE form of the system $k$:

$\left\\{\begin{array}[]{lcl}\dot{x}_{k}(t)&=&f_{k}(t,x_{k}(t),u_{k}(t))\\\
y_{k}(t)&=&g_{k}(t,x_{k}(t),u_{k}(t))\end{array}\right.$ (1)

Let $n_{in,tot}$ and $n_{out,tot}$ respectively be the total amount of inputs
$\sum_{k=1}^{n_{sys}}n_{in,k}$ and the total amount of outputs
$\sum_{k=1}^{n_{sys}}n_{out,k}$.

The total inputs and the total outputs vectors are simply concatenations of
input and output vectors of every system:

$\begin{array}[]{lclcl}\underline{\smash{u}}(t)&=&(u_{1}(t)^{T},\cdots,u_{n_{sys}}(t)^{T})^{T}&\in&L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{in,tot}})\\\
\underline{\smash{y}}(t)&=&(y_{1}(t)^{T},\cdots,y_{n_{sys}}(t)^{T})^{T}&\in&L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{out,tot}})\\\
\end{array}$ (2)

Finally, a tilde symbol $\tilde{}$ will be added to a functional quantity to
represent an element of its codomain. exempli gratia ,
$\underline{\smash{y}}\in L([t^{[N]},t^{[N+1]}[,\mathbb{R})$, so we can use
$\tilde{\underline{\smash{y}}}$ to represent an element of
$\mathbb{R}^{n_{out,tot}}$.

###  2.3 Extractors and rearrangement

In order to easily switch from global to local inputs, extractors are defined.
For $k\in[\\![1,n_{sys}]\\!]$, the extractor $E_{k}^{u}$ is the matrix defined
by (3).

$\begin{array}[]{rcccccl}E^{u}_{k}=\bigg{(}&0&\Big{|}&\Big{(}I_{n_{in,k}}\Big{)}&\Big{|}&0&\bigg{)}\\\
&\underbrace{\hskip 56.9055pt}&&\underbrace{\hskip
56.9055pt}&&\underbrace{\hskip 56.9055pt}&\\\ \lx@intercol
n_{in,k}\times\sum_{l=1}^{k-1}n_{in,l}\hfil\lx@intercol&&n_{in,k}\times
n_{in,k}&&\lx@intercol\hfil
n_{in,k}\times\sum_{l=k+1}^{n_{sys}}n_{in,l}\lx@intercol\\\ \end{array}$ (3)

where $\forall n\in\mathbb{N},\ I_{n}$ denotes the identity matrix of size $n$
by $n$.

The extractors enable to extract the inputs of a given system from the global
inputs vector with a relation of the form
$\tilde{u}_{k}=E^{u}_{k}\tilde{\underline{\smash{u}}}$. We have: $\forall
k\in[\\![1,n_{sys}]\\!],\ E^{u}_{k}\in M_{n_{in,k},n_{in,tot}}(\\{0,1\\})$.

A rearrangement operator will also be needed to handle concatenations of
outputs and output derivatives. For this purpose, we will use the
rearrangement matrix $R^{\underline{\smash{y}}}\in
M_{n_{out,tot},n_{out,tot}}(\\{0,1\\})$ defined blockwise in (4).

$\begin{array}[]{c}R^{\underline{\smash{y}}}=\left(R^{\underline{\smash{y}}}_{K,L}\right)_{\begin{subarray}{c}K\in[\\![1,\
2\ n_{sys}]\\!]\\\ L\in[\\![1,\ 2\ n_{sys}]\\!]\end{subarray}}\\\ \\\
\text{where}\\\ \\\
R^{\underline{\smash{y}}}_{K,L}=\left\\{\begin{array}[]{ll}I_{n_{out,K}}&\text{if}\
K\leqslant n_{sys}\ \text{and}\ L=2K-1\\\ I_{n_{out,K-n_{sys}}}&\text{if}\
K>n_{sys}\ \text{and}\ L=2(K-n_{sys})\\\ 0&\text{otherwise}\\\
\end{array}\right.\end{array}$ (4)

The $R^{\underline{\smash{y}}}$ operator makes it possible to rearrange the
outputs and output derivatives with a relation of the following form.

$\begin{array}[]{ccc}\left(\begin{array}[]{c}\\\
\tilde{\underline{\smash{y}}}\\\ \\\ \hline\cr\\\
\tilde{\dot{\underline{\smash{y}}}}\\\ \\\
\end{array}\right)=\left(\begin{array}[]{c}\tilde{y}_{1}\\\ \tilde{y}_{2}\\\
\vdots\\\ \tilde{y}_{n_{sys}}\\\ \hline\cr\tilde{\dot{y}}_{1}\\\
\tilde{\dot{y}}_{2}\\\ \vdots\\\ \tilde{\dot{y}}_{n_{sys}}\\\
\end{array}\right)=&\underbrace{\left(\begin{array}[]{ccccccc}\\!I_{n_{out,1}}\\!\\!&0&0&0&\cdots&0&0\\\
0&0&\\!\\!I_{n_{out,2}}\\!\\!&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&I_{n_{out,n_{sys}}}\\!\\!\\!&0\\\ \hline\cr
0&\\!I_{n_{out,1}}\\!\\!&0&0&\cdots&0&0\\\
0&0&0&\\!\\!I_{n_{out,2}}\\!&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&0&\\!\\!\\!I_{n_{out,n_{sys}}}\\!\\\
\end{array}\right)}&\left(\begin{array}[]{c}\tilde{y}_{1}\\\
\tilde{\dot{y}}_{1}\\\ \tilde{y}_{2}\\\ \tilde{\dot{y}}_{2}\\\ \vdots\\\
\tilde{y}_{n_{sys}}\\\ \tilde{\dot{y}}_{n_{sys}}\\\ \end{array}\right)\\\
&R^{\underline{\smash{y}}}&\\\ \end{array}$ (5)

###  2.4 Time discretization

In the context of co-simulation, the $g_{k}$ and $f_{k}$ functions in (1) are
usually not available directly. Thus, several co-simulation steps, so-called
”macro-steps”, are made between $t^{\text{init}}$ and $t^{\text{end}}$. Let’s
introduce the notations of the discrete version of the quantities introduced
in 2.2.

A macro step will be defined by its starting and ending times, respectively
denoted as $[t^{[N]},t^{[N+1]}]$ for the $N$th macro-step. The macro-steps
define a partition of the time-domain.

$\left\\{\begin{array}[]{lcl}[t^{\text{init}},t^{\text{end}}[&=&\displaystyle{\bigcup_{N=0}^{N_{\max}-1}}[t^{[N]},t^{[N+1]}[\\\
\\\ t^{[0]}&=&t^{\text{init}}\\\ t^{[N_{\max}]}&=&t^{\text{end}}\\\
\lx@intercol\hfil\forall N\in[\\![0,N_{\max}-1]\\!],\
<EMAIL_ADDRESS>(6)

Let $\delta t^{[N]}$ denote the size of the $N$th macro-step:

$\left\\{\begin{array}[]{l}\forall N\in[\\![0,N_{\max}-1]\\!],\ \delta
t^{[N]}=t^{[N+1]}-t^{[N]}>0\\\ \\\
\displaystyle{\sum_{N=0}^{N_{\max}-1}}\delta
t^{[N]}=t^{\text{end}}-t^{\text{init}}\end{array}\right.$ (7)

Let $\mathbb{T}$ denote the set of possible macro-steps.

$\mathbb{T}\overset{\Delta}{=}\\{[a,b[\ |\ t^{\text{init}}\leqslant
a<b\leqslant t^{\text{end}}\\}$ (8)

On a given macro-step $[t^{[N]},t^{[N+1]}[$, $N\in[0,N_{\max}]$, for all
systems, the restrictions of the piecewise equivalents of $u_{k}$ and $y_{k}$
will be denoted by $u_{k}^{[N]}$ and $y_{k}^{[N]}$ respectively. In case
several iterations are made on the same step, we will refer to the functions
by a left superscript index $m$. Finally, we will denote the coordinate of
these vectors with an extra subscript index.

$\begin{array}[]{llll}\lx@intercol\forall k\in[\\![1,n_{sys}]\\!],\ \forall
N\in[\\![0,N_{\max}]\\!],\ \forall m\in[0,m_{\max}(N)],\hfil\lx@intercol\\\
&\forall j\in[\\![1,n_{in,k}]\\!],&{}^{[m]}u_{k,j}^{[N]}&\in
L([t^{[N]},t^{[N+1]}[,\mathbb{R})\\\ &\forall
i\in[\\![1,n_{out,k}]\\!],&{}^{[m]}y_{k,i}^{[N]}&\in
L([t^{[N]},t^{[N+1]}[,\mathbb{R})\\\ \end{array}$ (9)

In (9), $m_{\max}(N)$ denotes the number of iterations (minus one) done on the
$N$th macro-step. $m_{max}(N)$ across $N$ can be plotted in order to see where
the method needed to proceed more or less iterations.

All derived notations introduced in this subsection can also be applied to the
total input and output vectors.

$\begin{array}[]{llll}\lx@intercol\forall N\in[\\![0,N_{\max}]\\!],\ \forall
m\in[\\![0,m_{\max}(N)]\\!],\hfil\lx@intercol\\\ &\forall\bar{\jmath}\in
n_{in,tot},&{}^{[m]}\underline{\smash{u}}_{\bar{\jmath}}^{[N]}&\in
L([t^{[N]},t^{[N+1]}[,\mathbb{R})\\\ &\forall\bar{\imath}\in
n_{out,tot},&{}^{[m]}\underline{\smash{y}}_{\bar{\imath}}^{[N]}&\in
L([t^{[N]},t^{[N+1]}[,\mathbb{R})\\\ \end{array}$ (10)

Indices $\bar{\imath}$ and $\bar{\jmath}$ in (10) will be called global
indices in opposition to the local indices $i$ and $j$ in (9).

###  2.5 Step function

Let $S_{k},\ k\in[\\![1,n_{sys}]\\!]$ be the ideal step function of the $k$th
system, that is to say the function which takes the system to its future state
one macro-step forward.

$S_{k}:\left\\{\begin{array}[]{lcl}\mathbb{T}\times
L([t^{\text{init}},t^{\text{end}}],\mathbb{R}^{n_{in,k}})\times\mathbb{R}^{n_{st,k}}&\rightarrow&\mathbb{R}^{n_{out,k}}\times\mathbb{R}^{n_{st,k}}\\\
(\tau,\ u_{k},\ \tilde{x})&\mapsto&S_{k}(\tau,\ u_{k},\
\tilde{x})\end{array}\right.$ (11)

In practice, the state vector $\tilde{x}$ will not be explicited. Indeed, it
will be embedded inside of the system $k$ and successive calls will either be
done:

  * $\bullet$

with $\tau$ beginning where the $\tau$ at the previous call of $S_{k}$ ended
(moving on),

  * $\bullet$

with $\tau$ beginning where the $\tau$ at the previous call of $S_{k}$ started
(step replay),

  * $\bullet$

with $\tau$ of the shape $[t^{\text{init}},t[$ with
$t\in]t^{\text{init}},t^{\text{end}}[$ (first step).

Moreover, the $u_{k}$ argument only needs to be defined on the domain $\tau$
(not necessary on $[t^{\text{init}},t^{\text{end}}[$). Thus, $S_{k}$ will not
be considered in the method, but the $\hat{S}_{k}$ function (practical step
function) defined hereafter will be considered instead. Despite $\hat{S}_{k}$
is not properly mathematically defined (the domain depends on the value of one
of the arguments: $\tau$ and some quantities are hidden: the states), it does
not lead to any problem, considering the hypotheses above.

$\begin{array}[]{c}\hat{S}_{k}:\left\\{\begin{array}[]{lcl}\mathbb{T}\times
L(\tau,\mathbb{R}^{n_{in,k}})&\mapsto&\mathbb{R}^{n_{out,k}}\\\
(\tau,u_{k})&\mapsto&\hat{S}_{k}(\tau,u_{k})\end{array}\right.\\\ \\\
\text{satisfying}\\\ \\\ \hat{S}_{k}([t^{[N]},t^{[N+1]}[,\
^{[m]}u_{k}^{[N]})=\ ^{[m]}y_{k}^{[N]}(t^{[N+1]})\end{array}$ (12)

The $\hat{S}_{k}$ function is the one available in practice, namely in the FMI
(Functionnal Mock-up Interface) standard.

###  2.6 Extended step function

The values of the output variables might not be sufficient for every co-
simulation scheme. It is namely the case for both classical IFOSMONDI and
IFOSMONDI-JFM. Indeed, the time-derivatives of the outputs are also needed.

Let $\hat{\hat{S}}_{k}$ be the extension of $\hat{S}_{k}$ returning both the
output values and derivatives.

$\begin{array}[]{c}\hat{\hat{S}}_{k}:\left\\{\begin{array}[]{lcl}\mathbb{T}\times
L(\tau,\mathbb{R}^{n_{in,k}})&\mapsto&\mathbb{R}^{n_{out,k}}\times\mathbb{R}^{n_{out,k}}\\\
(\tau,u_{k})&\mapsto&\hat{\hat{S}}_{k}(\tau,u_{k})\end{array}\right.\\\ \\\
\text{satisfying}\\\ \\\ \hat{\hat{S}}_{k}([t^{[N]},t^{[N+1]}[,\
^{[m]}u_{k}^{[N]})=\left(\begin{array}[]{c}{}^{[m]}y_{k}^{[N]}(t^{[N+1]})\\\
\\\ \displaystyle{\frac{d\
^{[m]}y_{k}^{[N]}}{dt}}(t^{[N+1]})\end{array}\right)\end{array}$ (13)

If the derivatives are not available in practice, a finite difference
approximation over $[t^{[N]},t^{[N+1]}[$ can be made (see
$\tilde{\hat{\hat{S}}}_{k}$ in [6]).

###  2.7 Connections

The connections between the systems will be denoted by a matrix filled with
zeros and ones, with $n_{out,tot}$ rows and $n_{in,tot}$ columns denoted by
$\Phi$. Please note that if each output is connected to exactely one input,
$\Phi$ is a square matrix. Moreover, it is a permutation matrix. Otherwise, if
an output is connected to several inputs, more than one $1$ appears at the
corresponding row of $\Phi$. In any case, there can neither be more nor less
than one $1$ on each column of $\Phi$ considering that an input can neither be
connected to none nor several outputs.

$\forall\bar{\imath}\in n_{out,tot},\ \forall\bar{\jmath}\in n_{in,tot},\
\Phi_{\bar{\imath},\bar{\jmath}}=\left\\{\begin{array}[]{ll}1&\text{if output
$\bar{\imath}$ is connected to input $\bar{\jmath}$}\\\
0&\text{otherwise}\end{array}\right.$ (14)

The dispatching will denote the stage where the inputs are generated from
their connected inputs, using the connections represented by $\Phi$.

$\tilde{\underline{\smash{u}}}=\Phi^{T}\tilde{\underline{\smash{y}}}$ (15)

The coupling function (16) will denote the absolute difference between
corresponding connected variables in a total input vector and a total output
vector. In other words, it represents the absolute error beween a total input
vector and the dispatching of a total output vector. The $\lambda$ subscript
does not correspond to a quantity, it is a simple notation inherited from a
Lagrange multipliers approach of systems coupling [13].

$g_{\lambda}:\left\\{\begin{array}[]{lcl}\mathbb{R}^{n_{in,tot}}\times\mathbb{R}^{n_{out,tot}}&\rightarrow&\mathbb{R}^{n_{in,tot}}\\\
(\tilde{\underline{\smash{u}}},\tilde{\underline{\smash{y}}})&\mapsto&|\tilde{\underline{\smash{u}}}-\Phi^{T}\tilde{\underline{\smash{y}}}|\end{array}\right.$
(16)

The coupling condition (17) is the situation where every output of the total
output vector corresponds to its connected input in the total input vector.

$g_{\lambda}(\tilde{\underline{\smash{u}}},\tilde{\underline{\smash{y}}})=0_{\mathbb{R}^{n_{in,tot}}}$
(17)

##  3 IFOSMONDI-JFM method

###  3.1 Modified extended step function

As in classical IFOSMONDI [6], the IFOSMONDI-JFM method preserves the $C^{1}$
smoothness of the interface variables at the communication times
$(t^{[N]})_{N\in[\\![1,N_{\max}-1]\\!]}$. Thus, when a time $t^{[N]}$ has been
reached, the input functions for every system will all satisfy the following
property:

$\begin{array}[]{llcl}\lx@intercol\forall k\in[\\![1,n_{sys}]\\!],\ \forall
m\in[\\![0,m_{\max}(N)]\\!],\hfil\lx@intercol\\\ \\\ \hskip
14.22636pt\hbox{\multirowsetup\bBigg@{4}\\{}&{}^{[m]}u_{k}^{[N]}(t^{[N]})&=&\
{}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}(t^{[N]})\\\ \\\ &\displaystyle{\frac{d\
^{[m]}u_{k}^{[N]}}{dt}}(t^{[N]})&=&\ \displaystyle{\frac{d\
^{[m_{\max}(N-1)]}u_{k}^{[N-1]}}{dt}}(t^{[N]})\\\ \end{array}$ (18)

The IFOSMONDI-JFM also represents the inputs as $3$rd order polynomial
(maximum) in order to satisfy the smoothness condition (18) and to respect
imposed values and derivatives at $t^{[N+1]}$ for every macro-step.

Knowing these constraints, it is possible to write a specification of the
practical step function $\hat{\hat{S}}_{k}$ in the IFOSMONDI-JFM case (also
applicable in the classical IFOSMONDI method):

$\zeta_{k}:\left\\{\begin{array}[]{lcl}\mathbb{T}\times\mathbb{R}^{n_{in,k}}\times\mathbb{R}^{n_{in,k}}&\mapsto&\mathbb{R}^{n_{out,k}}\times\mathbb{R}^{n_{out,k}}\\\
(\tau,\tilde{u}_{k},\tilde{\dot{u}}_{k})&\mapsto&\zeta_{k}(\tau,\tilde{u}_{k},\tilde{\dot{u}}_{k})\end{array}\right.\\\
$ (19)

where the three cases discussed in 2.5 have to be considered.

####  3.1.1 Case $1$: Moving on:

In this case, the last call to $\zeta_{k}$ was done with a $\tau\in\mathbb{T}$
ending at current $t^{[N]}$. In other words, the system $k$ ”reached” time
$t^{[N]}$. The inputs were, at this last call:
${}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}$.

To reproduce a behavior analog to the classical IFOSMONDI method, the inputs
${}^{[0]}u_{k}^{[N]}$ will be defined as the $2$nd order polynomial (or less)
satisfying the three following constraints:

$\begin{array}[]{lcl}{}^{[0]}u_{k}^{[N]}(t^{[N]})&=&{}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}(t^{[N]})\\\
\\\ \displaystyle{\frac{d\
^{[0]}u_{k}^{[N]}}{dt}}(t^{[N]})&=&\displaystyle{\frac{d\
^{[m_{\max}(N-1)]}u_{k}^{[N-1]}}{dt}}(t^{[N]})\\\ \\\
{}^{[0]}u_{k}^{[N]}(t^{[N+1]})&=&{}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}(t^{[N]})\\\
\end{array}$ (20)

The two first constraints guarantee the smoothness property (18), and the
third one minimizes the risk of out-of-range values (as in classical IFOSMONDI
method).

In this case, $\zeta_{k}$ in (19) is defined by the specification (21).

$\begin{array}[]{r}\zeta_{k}([t^{[N]},t^{[N+1]}[,\cdot,\cdot)=\hat{\hat{S}}_{k}([t^{[N]},t^{[N+1]}[,\
\underbrace{{}^{[0]}u_{k}^{[N]}})\\\ \text{computed with
\eqref{eq:inputs_def_case1_moving_on}}\end{array}$ (21)

$2$nd and $3$rd arguments of $\zeta_{k}$ are unused.

####  3.1.2 Case $2$: Step replay:

In this case, the last call to $\zeta_{k}$ was done with a $\tau\in\mathbb{T}$
starting at current $t^{[N]}$. In other words, the system did not manage to
reach the ending time of the previous $\tau$ (either because the method did
not converge, or because the step has been rejected, or another reason).

Two particular subcases have to be considered here: either the step we are
computing is following the previous one in the iterative method detailed after
this section, or the previous iteration has been rejected and we are trying to
re-integrate the step starting from $\tau$ with a smaller size $\delta
t^{[N]}$.

##### Subcase $2$.$1$: Following a previous classical step:

In this subcase, the last call of $\zeta_{k}$ was not only done with the same
starting time, but also with the same step ending time $t^{[N+1]}$. The inputs
were, at this last call: ${}^{[m-1]}u_{k}^{[N]}$ with $m\geqslant 1$, and
satisfied the two conditions at $t^{[N]}$ of (21).

The jacobian-free iterative method will ask for given input values
$\tilde{u}_{k}$ and time-derivatives $\tilde{\dot{u}}_{k}$ that will be used
as constraints at $t^{[N+1]}$, thus ${}^{[m]}u_{k}^{[N]}$ will be defined as
the $3$rd order polynomial (or less) satisfying the four following
constraints:

$\begin{array}[]{lclcl}{}^{[m]}u_{k}^{[N]}(t^{[N]})&=&{}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}(t^{[N]})&=&{}^{[m-1]}u_{k}^{[N]}(t^{[N]})\\\
\\\ \displaystyle{\frac{d\
^{[m]}u_{k}^{[N]}}{dt}}(t^{[N]})&=&\displaystyle{\frac{d\
^{[m_{\max}(N-1)]}u_{k}^{[N-1]}}{dt}}(t^{[N]})&=&\displaystyle{\frac{d\
^{[m-1]}u_{k}^{[N]}}{dt}}(t^{[N]})\\\ \\\
{}^{[m]}u_{k}^{[N]}(t^{[N+1]})&=&\tilde{u}_{k}\\\ \\\ \displaystyle{\frac{d\
^{[m]}u_{k}^{[N]}}{dt}}(t^{[N+1]})&=&\tilde{\dot{u}}_{k}\\\ \end{array}$ (22)

The two firsts constraints ensure the (18) smoothness property, and the third
and fourth ones will enable the iterative method to find the best values and
derivatives to satisfy the coupling condition.

In this subcase, $\zeta_{k}$ in (19) is defined by the specification (23).

$\begin{array}[]{r}\zeta_{k}([t^{[N]},t^{[N+1]}[,\tilde{u}_{k},\tilde{\dot{u}}_{k})=\hat{\hat{S}}_{k}([t^{[N]},t^{[N+1]}[,\
\underbrace{{}^{[m]}u_{k}^{[N]}})\\\ \text{computed with
\eqref{eq:inputs_def_case2_step_replay_subcase1}}\end{array}$ (23)

##### Subcase $2$.$2$: Re-integrate a step starting from $t^{[N]}$ but with
different $\delta t^{[N]}$ than at the previous call of $\zeta_{k}$:

In this subcase, current $t^{[N+1]}$ is different from $\sup{(\tau)}$ with
$\tau$ being the one used at the last call of $\zeta_{k}$.

As it shows that a step rejection just occured, we will simply do the same
than in case $1$, as if we were moving on from $t^{[N]}$. In other words, all
calls to $\zeta_{k}$ with $\tau$ starting at $t^{[N]}$ are ”forgotten”.

Please note that ${}^{[m_{\max}(N-1)]}u_{k}^{[N-1]}(t^{[N]})$ and
$\displaystyle{\frac{d\ ^{[m_{\max}(N-1)]}u_{k}^{[N-1]}}{dt}}(t^{[N]})$ can be
retreived using the values and derivatives constraints at $t^{[N]}$ of the
inputs at the last call of $\zeta_{k}$ thanks to the smoothness constraint
(18).

####  3.1.3 Case $3$: First step:

In this particular case, we will do the same as in the other cases, except
that we won’t impose any constraint for the time-derivative at
$t^{\text{init}}$. That is to say:

  * $\bullet$

at the first call of $\zeta_{k}$, we have $N=m=0$, we will only impose
${}^{[0]}u_{k}^{[0]}(t^{\text{init}})=\
^{[0]}u_{k}^{[0]}(t^{[1]})=u_{k}^{\text{init}}$ to have a zero order
polynomial satisfying the initial conditions $u_{k}^{\text{init}}$ (supposed
given),

  * $\bullet$

at the other calls, case $2$ will be used without considering the constraints
for the derivatives at $t^{\text{init}}$ (this will lower the polynomial’s
degrees). For (22), the first condition becomes
${}^{[m]}u_{k}^{[N]}(t^{\text{init}})=u_{k}^{\text{init}}$, the second one
vanishes, and the third ans fourth ones remain unchanged. For the subcase
$2$.$2$, it can be considered that
${}^{[m_{\max}(-1)]}u_{k}^{[-1]}(t^{\text{init}})=u_{k}^{\text{init}}$, and
$\displaystyle{\frac{d\ ^{[m_{\max}(-1)]}u_{k}^{[-1]}}{dt}}(t^{\text{init}})$
will not be needed as it is a time-derivative in $t^{\text{init}}$.

Finally, we have $\zeta_{k}$ defined in every case, wrapping polynomial inputs
computations and the integration done with $\hat{\hat{S}}_{k}$.

###  3.2 Iterative method’s callback function

The aim is to solve the co-simulation problem by using a jacobian-free version
of an iterative method that usually requires a jacobian computation (see 2.1).
Modern matrix-free versions of such algorithms make it possible to avoid
perturbating the systems and re-integrating them for every input, as done in
[13], in order to compute a finite-differences jacobian matrix. This saves a
lot of integrations over each macro-step and gains time.

Nevertheless, on every considered macro-step $\tau$, a function to be brought
to zero has to be defined. This so-called JFM’s callback (standing for
Jacobian-Free Method’s callback) presented hereafter will be denoted by
$\gamma_{\tau}$. In zero-order hold co-simulation, this function if often
$\tilde{\underline{\smash{u}}}-\Phi^{T}\tilde{\underline{\smash{y}}}$ (or
equivalent) where $\tilde{\underline{\smash{y}}}$ are the output at
$t^{[N+1]}$ generated by constant inputs $\tilde{\underline{\smash{u}}}$ over
$[t^{[N]},t^{[N+1]}[$.

In IFOSMONDI-JFM, we will only enable to change the inputs at $t^{[N+1]}$, the
smoothness condition at $t^{[N]}$ guaranteeing that the coupling condition
(17) remains satisfied at $t^{[N]}$ if it was satisfied before moving on to
the step $[t^{[N]},t^{[N+1]}[$. The time-derivatives will also be considered
in order to maintain $C^{1}$ smoothness, so the coupling condition (17) will
also be applied to these time-derivatives.

Finally, the formulation of the JFM’s callback for IFOSMONDI-JFM is:

$\gamma_{\tau}:\left\\{\begin{array}[]{lcl}\mathbb{R}^{n_{in,tot}}\times\mathbb{R}^{n_{in,tot}}&\rightarrow&\mathbb{R}^{n_{in,tot}}\times\mathbb{R}^{n_{in,tot}}\\\
\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)&\mapsto&\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)-\left(\begin{array}[]{cc}\Phi^{T}&0\\\
0&\Phi^{T}\end{array}\right)R^{\underline{\smash{y}}}\left(\begin{array}[]{c}\zeta_{1}\left(\tau,E^{u}_{1}\tilde{\underline{\smash{u}}},E^{u}_{1}\tilde{\dot{\underline{\smash{u}}}}\right)\\\
\vdots\\\
\zeta_{n_{sys}}\left(\tau,E^{u}_{n_{sys}}\tilde{\underline{\smash{u}}},E^{u}_{n_{sys}}\tilde{\dot{\underline{\smash{u}}}}\right)\\\
\end{array}\right)\end{array}\right.$ (24)

####  3.2.1 Link with the fixed-point implementation:

The formulation (24) can be used to represent the expression of the fixed-
point $\Psi_{\tau}$ function. The latter has been introduced in classical
IFOSMONDI algorithm [6] where a fixed-point method was used instead of a JFM
one. We can now rewrite a proper expression of $\Psi_{\tau}$ including the
time-derivatives.

$\Psi_{\tau}:\left\\{\begin{array}[]{lcl}\mathbb{R}^{n_{in,tot}}\times\mathbb{R}^{n_{in,tot}}&\rightarrow&\mathbb{R}^{n_{in,tot}}\times\mathbb{R}^{n_{in,tot}}\\\
\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)&\mapsto&\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)-\gamma_{\tau}(\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right))\\\
&&=\left(\begin{array}[]{cc}\Phi^{T}&0\\\
0&\Phi^{T}\end{array}\right)R^{\underline{\smash{y}}}\left(\begin{array}[]{c}\zeta_{1}\left(\tau,E^{u}_{1}\tilde{\underline{\smash{u}}},E^{u}_{1}\tilde{\dot{\underline{\smash{u}}}}\right)\\\
\vdots\\\
\zeta_{n_{sys}}\left(\tau,E^{u}_{n_{sys}}\tilde{\underline{\smash{u}}},E^{u}_{n_{sys}}\tilde{\dot{\underline{\smash{u}}}}\right)\\\
\end{array}\right)\end{array}\right.$ (25)

$\Psi_{\tau}$ was refered as $\Psi$ in [6] and did not include the derivatives
in its formulation, yet the smoothness enhancement done by the Hermite
interpolation led to an underlying use of these derivatives.

When the result of the $m$th iteration is available, a fixed-point iteration
on macro-step $\tau=[t^{[N]},t^{[N+1]}[$ is thus simply done by:

$\left(\begin{array}[]{c}{}^{[m+1]}\tilde{\underline{\smash{u}}}\\\
{}^{[m+1]}\tilde{\dot{\underline{\smash{u}}}}\end{array}\right):=\Psi_{\tau}(\left(\begin{array}[]{c}{}^{[m]}\tilde{\underline{\smash{u}}}\\\
{}^{[m]}\tilde{\dot{\underline{\smash{u}}}}\end{array}\right))$ (26)

###  3.3 First and last integrations of a step

The first iteration of a given macro-step $\tau\in\mathbb{T}$ is a particular
case to be taken into account. Considering the breakdown presented in
subsection 2.5, this corresponds to case $1$, case $2$ subcase $2$.$2$, case
$3$ first bullet point, and case $3$ second bullet point when falling into
subcase $2$.$2$.

All these cases have something in common: they denote calls to $\zeta_{k}$
using a $\tau$ argument that has never been used in a previous call of
$\zeta_{k}$. In these cases, the latter function is defined by (21).

For this reason, the first call of $\gamma_{\tau}$ for a given macro-step
$\tau$ will be done before applying the JFM. Then, every time the JFM will
call $\gamma_{\tau}$, the $(\zeta_{k})_{k\in[\\![1,n_{sys}]\\!]}$ functions
called by $\gamma_{\tau}$ will behave the same way.

Once the JFM method ends, if it converged, a last call to $\gamma_{\tau}$ is
made with the solution
$\big{(}(^{[m_{\max}(N)]}\tilde{\underline{\smash{u}}}^{[N]})^{T},\
(^{[m_{\max}(N)]}\tilde{\dot{\underline{\smash{u}}}}^{[N]})^{T}\big{)}^{T}$
for the systems to be in a good state for the next step (as explained in
subsection 2.5, the state of a system is hidden but affected at each call to a
step function).

###  3.4 Step size control

The step size control is defined with the same rule-of-thumbs than the one
used in [6]. The adaptation is not done on an error-based criterion such as in
[12], but with a predefined rule based on the convergence of the iterative
method (yes/no).

A reference step size $\delta t_{\text{ref}}\in\mathbb{R_{*}^{+}}$ is defined
for any simulation with IFOSMONDI-JFM method. It will either act as initial
macro-step size, and maximum step size. At some points, the method will be
allowed to reduce this step in order to help the convergence of the JFM.

The convergence criterion for the iterative method is defined by the rule
(27).

$\begin{array}[]{l}\text{Given}\
(\varepsilon_{\text{abs}},\varepsilon_{\text{rel}})\in(\mathbb{R_{+}^{*}})^{2},\\\
\text{convergence is reached when}\
\left|\gamma_{\tau}\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)\right|<\left|\left(\begin{array}[]{c}\tilde{\underline{\smash{u}}}\\\
\tilde{\dot{\underline{\smash{u}}}}\end{array}\right)\right|\varepsilon_{\text{rel}}+\left|\left(\begin{subarray}{c}1\\\
\vdots\\\ 1\end{subarray}\right)\right|\varepsilon_{\text{abs}}\end{array}$
(27)

When the iterative method does not converge on the step $[t^{[N]},t^{[N+1]}[$,
either because a maximum number of iterations is reached or for any other
reason (linear search does not converge, a Krylov internal method finds a
singular matrix, …), the step will be rejected and retried on the half (28).
Otherwise, once the method converged on $[t^{[N]},t^{[N+1]}[$, the next
integration step $\tau$ tries to increase the size of $30\%$, without
exceeding $\delta t_{\text{ref}}$.

Once the iterative method exits on $\tau_{\text{old}}$, the next step
$\tau_{\text{new}}$ is defined by:

$\small\tau_{\text{new}}=\left\\{\begin{array}[]{lr}\lx@intercol\left[\sup(\tau_{\text{old}}),\min\\!\bigg{\\{}t^{\text{end}},\sup(\tau_{\text{old}})+\max\\!\Big{\\{}\delta
t_{\text{ref}},1.3\
\big{(}\sup(\tau_{\text{old}})-\inf(\tau_{\text{old}})\big{)}\Big{\\}}\bigg{\\}}\right[\hfil\lx@intercol\\\
&\text{if convergence \eqref{eq:convergence_criterion} was reached}\\\ \\\
\lx@intercol\left[\inf(\tau_{\text{old}}),\inf(\tau_{\text{old}})+\displaystyle{\frac{\sup(\tau_{\text{old}})-\inf(\tau_{\text{old}})}{2}}\right[\hfil\lx@intercol\\\
&\text{otherwise (divergence)}\\\ \end{array}\right.$ (28)

When $\varepsilon_{\text{abs}}=\varepsilon_{\text{rel}}$, these values will be
denoted by $\varepsilon$.

##  4 Note on the implementation

Our implementation is based on an orchestrator-worker architecture, where
$n_{sys}+1$ processes are involved. One of them is dedicated to global
manipulations: the orchestrator. It is not responsible of any system and only
deals with global quantities (such as the time, the step $\tau$, the
$\tilde{\underline{\smash{u}}}$ and $\tilde{\underline{\smash{y}}}$ vectors
and the corresponding time-derivatives, and so on). The $n_{sys}$ remaining
processes, the workers, are reponsible of one system each. They only deal with
local quantities related to the system they are responsible of.

###  4.1 Parallel evaluation of $\gamma_{\tau}$ using MPI

An evaluation of $\gamma_{\tau}$ consists in evaluations of the $n_{sys}$
functions $(\zeta_{k})_{k\in[\\![1,n_{sys}]\\!]}$, plus some manipulations of
vectors and matrices (24). An evaluation of a single $\zeta_{k}$ for a given
$k\in[\\![1,n_{sys}]\\!]$ consists in polynomial computations and an
integration (21) (23) through a call of the corresponding $\hat{\hat{S}}_{k}$
function (13).

A single call to $\gamma_{\tau}$ can be evaluated parallely by $n_{sys}$
processes, each of them carying out the integration of one of the systems. To
achieve this, the MPI standard (standing for Message Passing Interface has
been used, as the latter provides routine to handle multi-process
communications of data.

As the $k$th system only needs $E^{u}_{k}\tilde{\underline{\smash{u}}}$ and
$E^{u}_{k}\tilde{\dot{\underline{\smash{u}}}}$ (see (3)) among
$\tilde{\underline{\smash{u}}}$ and $\tilde{\dot{\underline{\smash{u}}}}$, the
data can be send in an optimized manner from an orchestrator process to
$n_{sys}$ workers by using the `MPI_Scatterv` routine.

Analogously, each worker process will have to communicate their contribution
both to the outputs and their derivatives (assembling the block vector at the
right of the expression (24)). This can be done by using the `MPI_Gatherv`
routine.

Finally, global quantities such as $\tau$, $m$, the notifications of statuses
and so on can be done easily thanks to the `MPI_Broadcast` routine.

###  4.2 Using PETSc for the JFM

PETSc [2, 3] is a library used for parallel numerical computations. In our
case, the several matrix-free versions of the Newton method and variants
implemented in PETSc were very attractive. Indeed, the flexibility of this
library at runtime enables the use of command-line arguments to control the
resolution: `-snes_mf` orders the use of a matrix-free non-linear solver,
`-snes_type newtonls`, `anderson` [1] and `ngmres` [10] are various usable
solving methods that can be used as JFMs, `-snes_atol`, `-snes_rtol` and
`-snes_max_it` control the convergence criterion, `-snes_converged_reason`,
`-snes_monitor` and `-log_view` produce information and statistics about hte
run, …

This subsection proposes a solution to use these PETSc implementations in a
manner that is compliant with the parallel evaluation of the JFM’s callback
(24). This implementation has been used to generate the results of section 5.

First of all, PETSc needs a view on the environment of the running code: the
processes, and their relationships. In our case, the $n_{sys}+1$ processes of
the orchestrator-worker architecture are not dedicated to the JFM. Thus, PETSc
runs on the orchestrator process only. In terms of code, this can be done by
creating PETSc objects referring to `PETSC_COMM_SELF` communicator on the
orchestrator process, and creating no PETSc object on the workers.

The callback $\gamma_{\tau}$ implements internally the communications with the
workers, and is given to the PETSc `SNES` object. The `SNES` non-linear solver
will call this callback blindly, and the workers will be triggered behind the
scene for integrations, preceded by the communications of the
$\big{(}(^{[m_{\max}(N)]}\tilde{\underline{\smash{u}}}^{[N]})^{T},\
(^{[m_{\max}(N)]}\tilde{\dot{\underline{\smash{u}}}}^{[N]})^{T}\big{)}^{T}$
values asked by the `SNES` and followed by the gathering of the outputs and
related derivatives. The latters are finally returned to PETSc by the callback
on the orchestrator side, after reordering and dispatching them as in (24).

###  4.3 JFM’s callback implementation

In this section, a suggestion of implementation is proposed for the
$\gamma_{\tau}$ function, both on the orchestrator side and on the workers
side. Precisions about variables in the snippets are given below them.

The aim is not to show the code that has been used to generate the results of
section 5, but to figure out how to combine the PETSc and MPI standard (PETSc
being based on MPI) to implement a parallel evaluation of $\gamma_{\tau}$.

By convention, the process of rank $0$ is the orchestrator, and any process of
rank $k\in[\\![1,n_{sys}]\\!]$ will be responsible of system $k$.

Snippet 1: JFM’s callback on the orchestrator side ($\gamma_{\tau}$)

⬇

1PetscErrorCode JFM_callback(SNES /* snes */, Vec u_and_du, Vec res, void
*ctx_as_void)

2{

3 MyCtxType *ctx = (MyCtxType*)ctx_as_void;

4 const int order = DO_A_STEP;

5 PetscScalar const * pscalar_u_and_du;

6 PetscScalar * pscalar_res;

7 size_t k;

8

9 // conversion PetscScalar -> C double

10 VecGetArrayRead(u_and_du, &pscalar_u_and_du);

11 for (k = 0; k < ctx->n_in_tot * 2; k++)

12 ctx->double_u_and_du[k] = (double)(pswork_x[k]);

13 VecRestoreArrayRead(u_and_du, &pscalar_u_and_du);

14

15 // Notify workers that we want them to run,

16 // and telling them what tau is

17 MPI_Bcast(&order, 1, MPI_INT, 0, MPI_COMM_WORLD);

18 MPI_Bcast(&(ctx->t_N), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

19 MPI_Bcast(&(ctx->t_Np1), 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

20

21 // Apply extrators and communicate at the same time:

22 // values, then derivatives

23 MPI_Scatterv(ctx->double_u_and_du, ctx->in_sizes, ctx->in_offsets,
MPI_DOUBLE, NULL, 0, MPI_DOUBLE, 0, MPI_COMM_WORLD);

24 MPI_Scatterv(ctx->double_u_and_du + ctx->n_in_tot, ctx->in_sizes,
ctx->in_offsets, MPI_DOUBLE, NULL, 0, MPI_DOUBLE, 0, MPI_COMM_WORLD);

25

26 /* Workers proceed integration here */

27

28 // Assemble vector R^{\bar{y}} * (\zeta_1^T, ... \zeta_2^T)^T directly

29 // while communicating values and derivatives

30 MPI_Gatherv_outputs(MPI_IN_PLACE, 0, MPI_DOUBLE, ctx->work1_n_out_tot,
ctx->out_sizes, ctx->out_offsets, MPI_DOUBLE, 0, MPI_COMM_WORLD);

31 MPI_Gatherv_outputs(MPI_IN_PLACE, 0, MPI_DOUBLE, ctx->work2_n_out_tot,
ctx->out_sizes, ctx->out_offsets, MPI_DOUBLE, 0, MPI_COMM_WORLD);

32

33 // Dispatching (equivalent of [[Phi^T, 0], [0, Phi^T]])

34 dispatch(ctx->work1_n_out, ctx->out_sizes, ctx->n_sys,

35 ctx->double_res, ctx->in_sizes, ctx->connections);

36 dispatch(ctx->work2_n_out, ctx->out_sizes, ctx->n_sys,

37 ctx->double_res + ctx->n_in_tot, ctx->in_sizes, ctx->connections);

38

39 // Difference between original entries and permuted outputs

40 for (k = 0; k < ctx->n_in_tot * 2; k++)

41 ctx->double_res[k] = ctx->double_u_and_du[k] - ctx->double_res;

42

43 // conversion C double -> PetscScalar

44 VecGetArray(res, &pscalar_res);

45 for (k = 0; k < ctx->n_in_tot * 2; k++)

46 pswork_f[k] = (PetscScalar)(ctx->double_res[k]);

47 VecRestoreArray(res, &pscalar_res);

48

49 return 0;

50}

In the code snippet 1, the function `JFM_callback` is the one that is given to
the PETSc `SNES` object with `SNESSetFunction`. The context pointer `ctx` can
be anything that can be used to have access to extra data inside of this
callback. The principle is: when `SNESSolve` is called, the callback function
which has been given to the `SNES` object will be called an unknown number of
times. For this example, we suggested a context structure `MyCtxType` at least
containing:

  * $\bullet$

`t_N`, `t_Np1` the boundary times of $\tau$, id est $t^{[N]}$ and $t^{[N+1]}$
(as `double` each),

  * $\bullet$

`n_in_tot` the total number of inputs $n_{in,tot}$ (as `size_t`),

  * $\bullet$

`double_u_and_du` an array dedicated to the storage of
$(\tilde{\underline{\smash{u}}}^{T},\tilde{\dot{\underline{\smash{u}}}}^{T})^{T}$
(as `double*`),

  * $\bullet$

`in_sizes` the array containing the number of inputs for each process
$(n_{in,k})_{k\in[\\![0,n_{sys}]\\!]}$ including process $0$ (with the
convention $n_{in,0}=0$) (as `int*`),

  * $\bullet$

`in_offsets` the memory displacements
$\left(\sum_{l=1}^{k}n_{in,l}\right)_{k\in[\\![0,n_{sys}]\\!]}$ for inputs
scattering for each process (as `int*`),

  * $\bullet$

`work1_n_out_tot` and `work2_n_out_tot` two arrays of size $n_{out,tot}$ for
temporary storage (as `double*`),

  * $\bullet$

`out_sizes` and `out_offsets` two arrays analogous to `in_sizes` and
`in_offsets` respectively, considering the outputs,

  * $\bullet$

`n_sys` tot number of systems $n_{sys}$ (as `size_t`),

  * $\bullet$

`double_res` an array of size $2\ n_{in,tot}$ dedicated to the storage of the
result of $\gamma_{\tau}$ (as `double*`), and

  * $\bullet$

`connections` any structure to represent the connections between the systems
$\Phi^{T}$ (a full matrix might be a bad idea as $\Phi$ is expected to be very
sparse).

Finally, `dispatch` is expected to be a function processing the dispatching
(15) of the values given in its first argument into the array pointed by its
fourth argument.

On the workers side, the running code section is the one in the snippet 2.

Snippet 2: JFM’s callback on the worker side ($\zeta_{k}$ and communications)

⬇

1/* ... */

2

3while (1)

4{

5 MPI_Bcast(&order, 1, MPI_INT, 0, me__->comm);

6 if (order != DO_A_STEP)

7 break;

8

9 // get tau

10 MPI_Bcast(&t_N, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

11 MPI_Bcast(&t_Np1, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

12

13 // receive relevant inputs and derivatives for this system

14 MPI_Scatterv(NULL, ctx->in_sizes, ctx->in_offsets, MPI_DOUBLE,

15 sys_inputs, sys_n_in, MPI_DOUBLE, 0, me__->comm);

16 MPI_Scatterv(NULL, ctx->in_sizes, ctx->in_offsets, MPI_DOUBLE,

17 sys_dinputs, sys_n_in, MPI_DOUBLE, 0, me__->comm);

18

19 /* integration: */

20 zeta_do_a_step(t_N, t_Np1, inputs, sys_dinputs, // [in]

21 sys_outputs, sys_doutputs); // [out]

22

23 // send the outputs and derivatives (results of zeta)

24 MPI_Gatherv_outputs(sys_outputs, sys_n_out, MPI_DOUBLE,

25 NULL, NULL, NULL, MPI_DOUBLE, 0, MPI_COMM_WORLD);

26 MPI_Gatherv_outputs(sys_doutputs, sys_n_out, MPI_DOUBLE,

27 NULL, NULL, NULL, MPI_DOUBLE, 0, MPI_COMM_WORLD);

28}

29

30/* ... */

Please note that the orchestrator process has to explicitely send an order
different from `DO_A_STEP` (with `MPI_Bcast`) to notify the workers that the
callback will not be called anymore.

Nonetheless, this order might not be send right after the call to `SNESSolve`
on the orchestrator side. Indeed, if the procedure converged, a last call has
to be made explicitely in the orchestrator (see 3.3).

An other explicit call to `JFM_callback` should also be explicitly made on the
orchestrator side before the call of `SNESSolve` (as also explained in 3.3).

Finally, figure 4.3 presents a schematic view of these two snippets running
parallely.

Figure 1: Workflow of the callback function called by SNESSolve: example with
$n_{sys}=2$ (external first call to the callback is supposed to be already
made before SNESSolve is called)

##  5 Results on a test cases

Difficulties may appear in a co-simulation problem when the coupling is not
straightforward. Some of the most difficult cases to solve are the algebraic
coupling (addressed in [8]) arising from causal conflicts, and the multiple
feed-through, id est the case where outputs of a system linearly depend on its
inputs, and the connected system(s) have the same behavior. In some case, this
may lead to a non-contractant $\Psi_{\tau}$ function. This section presents a
test case we designed, belonging to this second category. The fixed-point
convergence can be accurately analyzed so that its limitations are
highlighted.

###  5.1 Test case presentation

Figure 2: Mass spring damper with damping reaction modelled with Simcenter
Amesim - Parameters are above, variables are below

Figure 5.1 represents a $1$-mass test case with a classical mechanical
coupling on force, velocity and position. These coupling quantities are
respectively denoted by $f_{c}$, $v_{c}$ and $x_{c}$. The component on the
right represents a damper with a massless plate, computing a velocity (and
integrating it to compute a displacement) by reaction to a force input.

We propose the parameters values in table 5.1.

Table 1: Parameters and initial values of the test case model Notation | Description | Value
---|---|---
$M_{L}$ | Mass of the body in $(S_{1})$ |  $1$ kg
$K_{SD}$ | Spring rate of the spring in $(S_{1})$ |  $1$ N/m
$D_{SD}$ | Damper rating of the damper in $(S_{1})$ |  $1$ N/(m/s)
$D_{D}$ | Damper rating of the damper in $(S_{2})$ | $\in[0.01,4]$
$x_{L}(0)$ | Initial position of the body in $(S_{1})$ |  $0$ m
$v_{L}(0)$ | Initial velocity of the body in $(S_{1})$ |  $0$ m/s
$x_{D}(0)$ | Initial position of the plate in $(S_{2})$ |  $0$ m
$t^{\text{init}}$ | Initial time |  $0$ s
$t^{\text{end}}$ | Final time |  $10$ s

All variables will be denoted by either $f$, $v$ or $x$ (corresponding to
forces, velocities and positions, respectively) with an index specifying its
role in the model (see figure 5.1).

The predefined force $f_{L}$ is a $C^{\infty}$ function starting from $5$ N
and definitely reaching $0$ N at $t=2$ s. The expression of $f_{L}$ is (29)
and the visualization of it is presented on figure 5.1.

Figure 3: Predefined force $f_{L}$

$f_{L}:\left\\{\begin{array}[]{lcl}[0,10]&\rightarrow&[0,5]\\\
t&\mapsto&\left\\{\begin{array}[]{ll}\displaystyle{\frac{5}{e^{-1}}}e^{\left(\left(\displaystyle{\frac{t}{2}}\right)^{2}-\
1\right)^{-1}}&\text{if}\ t<2\\\ \\\ 0&\text{if}\ t\geqslant 2\\\
\end{array}\right.\end{array}\right.$ (29)

The expected behavior of the model is presented in table 5.1 referring to
conventionnal directions of figure 5.1.

Figure 4: Test model visualized with Simcenter Amesim

Table 2: Main stages of a simulation of the test case model Stage | Body | Plate | Description
---|---|---|---
displacement | displacement
$1$ | front | front | Positive $f_{L}$ pushes everything
$2$ | back | front |  | The spring pushes the body
---
backward as it is close to the plate
$3$ | back | back |  | The spring pulls the plate backwards as
---
the body is moving backward with inertia
$4$ | front | back |  | The spring pulls the body forward as the inertia
---
made it go too far in the backward direction
$5$ | front | front |  | The body is still moving frontward with inertia,
---
so the compressed spring pushed the plate forward

The behavior presented in table 5.1 might slightly change while parameter
$D_{D}$ changes (all other parameters being fixed, see table 5.1).

###  5.2 Equations and eigenvalues of the fixed-point callback $\Psi_{\tau}$

Second Newton’s law gives:

$\begin{array}[]{lcl}\dot{v}_{L}&=&(f_{L}+f_{SD})M_{L}^{-1}\\\
\dot{x}_{L}&=&v_{L}\end{array}$ (30)

and the spring and damper forces can be expressed the following way:

$\begin{array}[]{lcl}f_{SD}&=&K_{SD}(x_{C}-x_{L})+D_{SD}(v_{C}-v_{L})\\\
f_{C}&=&-f_{SD}\\\ f_{D}&=&-D_{D}(0-v_{C})\\\ f_{C}&=&f_{D}\\\
v_{C}&=&\nicefrac{{f_{C}}}{{D_{D}}}\end{array}$ (31)

leading to the following expressions of the coupled systems:

$\begin{array}[]{l}(S_{1}):\left\\{\begin{array}[]{ccl}\left(\begin{array}[]{c}\dot{v}_{L}\\\
\dot{x}_{L}\end{array}\right)&=&\left(\begin{array}[]{cc}\frac{-D_{SD}}{M_{L}}&\frac{-K_{SD}}{M_{L}}\\\
1&0\end{array}\right)\left(\begin{array}[]{c}v_{L}\\\
x_{L}\end{array}\right)+\left(\begin{array}[]{cc}\frac{D_{SD}}{M_{L}}&\frac{K_{SD}}{M_{L}}\\\
0&0\end{array}\right)\left(\begin{array}[]{c}v_{C}\\\
x_{C}\end{array}\right)+\left(\begin{array}[]{c}\frac{f_{L}}{M_{L}}\\\
0\end{array}\right)\\\ f_{C}&=&\left(D_{SD}\ \
K_{SD}\right)\left(\begin{array}[]{cc}v_{L}\\\
x_{L}\end{array}\right)+\left(-D_{SD}\ \
-K_{SD}\right)\left(\begin{array}[]{cc}v_{C}\\\
x_{C}\end{array}\right)\end{array}\right.\\\ \\\
(S_{2}):\left\\{\begin{array}[]{ccl}\dot{x}_{D}&=&0\ x_{D}+\frac{1}{D_{D}}\
f_{C}\\\ \left(\begin{array}[]{c}v_{C}\\\
x_{C}\end{array}\right)&=&\left(\begin{array}[]{c}0\\\
1\end{array}\right)x_{D}+\left(\begin{array}[]{c}\frac{1}{D_{D}}\\\
0\end{array}\right)f_{C}\end{array}\right.\end{array}$ (32)

At a given time $t$, we can state the jacobian of $\Psi_{\tau}$ introduced in
(25) using the expressions of the coupling quantities (32). Indeed, the output
variables got at a call are at the same time than the one at which the imposed
inputs are reached (end of the macro-step) thanks to the definitions of
$\zeta_{k}$.

$J_{\Psi_{\tau}}(\left(\begin{array}[]{c}f_{C}\\\ v_{C}\\\ x_{c}\\\
\dot{f_{C}}\\\ \dot{v_{C}}\\\
\dot{x_{c}}\end{array}\right))=\left(\begin{array}[]{ccccccccc}\hline\cr\vrule\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&&-D_{SD}&-K_{SD}&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&&0&0\\\
\hline\cr\nicefrac{{1}}{{D_{D}}}&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&&0&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil 0\hfil\lx@intercol\vrule\lx@intercol\\\
0&&\vrule\lx@intercol\hfil 0\hfil\lx@intercol&\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&&0&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil 0\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\vspace{-3mm}\hfil\\\
\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol\vrule\lx@intercol&&-D_{SD}&-K_{SD}\\\
\hline\cr\vrule\lx@intercol\hfil\text{Block}\hfil\lx@intercol\vrule\lx@intercol&&\nicefrac{{1}}{{D_{D}}}&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil 0\hfil\lx@intercol\vrule\lx@intercol\\\
\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&&0&&\vrule\lx@intercol\hfil
0\hfil\lx@intercol&\lx@intercol\hfil 0\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\end{array}\right)$ (33)

The framed zeros are ”by-design” zeros: indeed, systems never produce outputs
depending on inputs given to other systems. The block called ”Block” in (33)
depends on the methods used to retrieve the time-derivatives of the coupling
quantities (see (13) and its finite differences version). Nevertheless, this
block does not change the eigenvalues of $J_{\Psi_{\tau}}$ as it is a block-
triangular matrix. Indeed, the characteristic polynomial of $I_{6}-\lambda
J_{\Psi_{\tau}}$ is the product of the determinant of the two $3\times 3$
blocks on the diagonal of $I_{6}-\lambda J_{\Psi_{\tau}}$. The eigenvalues of
$J_{\Psi}$ are:

$0,\ +1\\!\textit{i}\sqrt{\frac{D_{SD}}{D_{D}}},\
-1\\!\textit{i}\sqrt{\frac{D_{SD}}{D_{D}}}\hskip 14.22636pt\text{(each with a
multiplicity of $2$)}$ (34)

Hence, the following relation between the parameters and the spectral radius
can be shown (given $D_{D}>0$ and $D_{SD}=1>0$):

$\varrho\left(J_{\Psi_{\tau}}\right)\left\\{\begin{array}[]{lcl}<1&\
\text{if}\ D_{SD}<D_{D}\\\ \geqslant 1&\ \text{if}\ D_{SD}\geqslant
D_{D}\end{array}\right.$ (35)

We can thus expect that the classical IFOSMONDI co-simulation algorithm based
on a fixed-point method [6] cannot converge on this model when the damping
ratio of the component on the right of the model (see figure 5.1) is smaller
than the damping ratio of the spring-damper component.

We will process several simulations with different values of $D_{D}$ leading
to different values of $\varrho(J_{\Psi_{\tau}})$. These values and the
expected movement of the body of the system is plotted in figure 5.2.

Figure 5: Displacement of the mass ($x_{L}$) for different damping ratios of
the right damper ($D_{D}$) simulated on a monolithic model (without co-
simulation). Associated spectral radii of $J_{\Psi}$ are recalled for futher
coupled formulations.

###  5.3 Results

As the PETSc library enables to easily change the parameters of the JFM (as
explained in subsection 4.2), three methods have been used in the simulations:

  * $\bullet$

NewtonLS: a Newton based non-linear solver that uses a line search,

  * $\bullet$

Ngmres: the non-linear generalized minimum residual method [10], and

  * $\bullet$

Anderson: the Anderson mixing method [1]

First of all, simulations have been processed with all these JFMs (with
parameters exhaustively defined in appendix A) within IFOSMONDI-JFM, the
classical IFOSMONDI algorithm (denoted hereafter as ”Fixed-point”), and the
original explicit zero-order hold co-simulation method (sometimes referred to
as non-iterative Jacobi). The error is defined as the mean of the normalized
$L^{2}$ errors on each state variable of both systems on the whole
$[t^{\text{init}},t^{\text{end}}]$ domain. The reference is the monolithic
simulation (of the non-coupled model) done with Simcenter Amesim. Such errors
are presented for a contractant case ($D_{D}=4$ N, so
$\varrho(J_{\Psi_{\tau}})=0.5$) in figure 5.3. For a non-contractant case
($D_{D}=0.64$ N, so $\varrho(J_{\Psi_{\tau}})=1.25$), analog plots are
presented in figure 5.3.

Figure 6: Error accross $\delta t_{\text{ref}}$ with different methods on a
contractant case ($D_{D}=4.0$, $\rho(J_{Psi})=0.5$) - NewtonLS, Ngmres and
Anderson are matrix-free iterative methods used with the IFOSMONDI-JFM
algorithm, Fixed-point is the classical IFOSMONDI algorithm, and Explicit ZOH
is the non-iterative zero-order hold fixed-step co-simulation Figure 7: Error
accross $\delta t_{\text{ref}}$ with different methods on a non-contractant
case ($D_{D}=0.64$, $\rho(J_{Psi})=1.25$) - NewtonLS, Ngmres and Anderson are
matrix-free iterative methods used with the IFOSMONDI-JFM algorithm

As expected, the simulations failed (diverged) with fixed-point method
(classical IFOSMONDI) for the non-contractant case. Moreover, the values given
to the system were too far from physically-possible values with the explicit
ZOH co-simulation algorithm, so the internal solvers of systems $(S_{1})$ and
$(S_{2})$ failed to integrate. These are the reason why these two methods lead
to no curve on figure 5.3.

Nonetheless, the three versions of IFOSMONDI-JFM algorithm keep producing
reliable results with an acceptable relative error (less than $1\%$) when
$\delta t_{\text{ref}}\geqslant 0.1$ s.

On figures 5.3 and 5.3, IFOSMONDI-JFM method seems to solve the problem with a
good accuracy regardless of the value of the damping ratio $D_{D}$. In order
to confirm that, several other values have been tried: the ones for which the
solution has been computed and plotted in figure 5.2. The error is presented,
but also the number of iterations and the number of integrations (calls to
$\zeta_{k}$, i.e. calls to $\gamma_{\tau}$ for IFOSMONDI-JFM or to
$\Psi_{\tau}$ for classical IFOSMONDI). Although for the fixed-point case
(classical IFOSMONDI) the number of iteration is the same than the number of
integration, for the IFOSMONDI-JFM algorithm the number of iterations is the
one of the underlying non-linear solver (NewtonLS, Ngmres or Anderson), and
there might be a lot more integrations than iterations of the non-linear
method. These results are presented in figure 5.3.

Figure 8: Total number of iterations, integrations, and error across spectral
radius of $J_{\Psi}$ for different methods (Fixed-point corresponds to
classical IFOSMONDI algorithms, and all other methods are used with the
IFOSMONDI-JFM version). All co-simulation ran with $\varepsilon=10^{-4}$ and
$\delta t_{\text{ref}}=10^{-2}$

As expected, the threshold of $\varrho(J_{\Psi_{\tau}})=1$ (id est
$D_{D}=D_{SD}=1$) is critical for the fixed-point method. The IFOSMONDI-JFM
method not only can overpass this threshold, but no significant extra
dificulty appears to solve the problem in the non-contractant cases, except
for the Ngmres non-linear solver (which failed to converge with $D_{D}=0.01$,
so with $\varrho(J_{\Psi_{\tau}})=10$). However, regarding the Ngmres method,
the variant that uses line search converges in all cases. Eventhough the
latter requires more integrations than other JFMs, it is more robust to high
values of $\varrho(J_{\Psi_{\tau}})$. The parameters of this line search are
detailed on table A in appendix A.

The NewtonLS and Anderson methods show a slightly bigger error on this
”extreme” case of $\varrho(J_{\Psi_{\tau}})=10$, yet it stays under $0.001\%$
which is completely acceptable.

Among those two JFMs (NewtonLS and Anderson), the trend that can be observed
on figure 5.3 shows that NewtonLS is always more accurate than Anderson, yet
it always requires a bigger amount of integrations. We can stand that
IFOSMONDI-JFM is more accuracy-oriented on this model when it is based on the
NewtonLS JFM, and more speed-oriented on this model when it is based on the
Anderson JFM (for the same $\delta t_{\text{ref}}$ and $\varepsilon$). For
high values of $\varrho(J_{\Psi_{\tau}})$, accuracy-oriented simulations are
achieved thanks to the Ngmres JFM with line search more than the NewtonLS one.

Finally, smaller errors are obtained with IFOSMONDI-JFM and with less
iterations than classical IFOSMONDI. Yet, the time consumption is directly
linked with the number of integrations, not with the number of iterations of
the underlying non-linear solver. The total number of integrations does not
increase across the problem difficulty (increasing with
$\varrho(J_{\Psi_{\tau}})$), and the non-linear methods within IFOSMONDI-JFM
do not even require more integrations that the fixed-point one for most of the
values of $D_{D}$ for which the classical IFOSMONDI algorithm does not fail.

##  6 Conclusion

IFOSMONDI-JFM method not only enables to solve problems that were impossible
to solve with the classical IFOSMONDI method, it also requires less iterations
to converge on the test case of section 5 when the parameterization enables
both methods to solve the problem.

Despite a number of integration greater than one for every iteration (contrary
to the classical IFOSMONDI algorithm), IFOSMONDI-JFM does not require a lot
more integrations than classical IFOSMONDI. In most of the cases, for the
considered test case, IFOSMONDI-JFM even requires less integrations than
classical IFOSMONDI, and the resulting solution is always more accurate (for
the same $\delta t_{\text{ref}}$ and $\varepsilon$). The matrix-free aspect of
the underlying solvers used with IFOSMONDI-JFM are one of the causes of the
small amount of integrations.

The IFOSMONDI-JFM algorithm takes advantages from the $C^{1}$ smoothness of
classical IFOSMONDI algorithm [6] without the delay it implies in [5] (thanks
to its iterative aspect), the coupling constraint is satisfied both at left
and right of every communication time thanks to the underlying non-linear
solvers of PETSc [2]. The iterative part does not need a finite differences
estimation of the jacobian matrix like in [13] or a reconstruction of it like
in [14].

The resulting algorithm even solves co-simulation problems for which the
fixed-point formulation would involve a non-contractant coupling function
$\Psi_{\tau}$.

Finally, the test case introduced in 5.1 can be reused to test the robustness
of various co-simulation methods as the model is relatively simple and the
difficulty can easily be increased or decreased in a quantifiable way.

## References

  * [1] Anderson, D.G.M.: Iterative Procedures for Nonlinear Integral Equations. J. ACM 12, 547–560 (1965)
  * [2] Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W.D., Karpeyev, D., Kaushik, D., Knepley, M.G., May, D.A., McInnes, L.C., Mills, R.T., Munson, T., Rupp, K., Sanan, P., Smith, B.F., Zampini, S., Zhang, H., Zhang, H.: PETSc Web page. https://www.mcs.anl.gov/petsc (2019), https://www.mcs.anl.gov/petsc
  * [3] Balay, S., Gropp, W.D., McInnes, L.C., Smith, B.F.: Efficient management of parallelism in object oriented numerical software libraries. In: Arge, E., Bruaset, A.M., Langtangen, H.P. (eds.) Modern Software Tools in Scientific Computing. pp. 163–202. Birkhäuser Press (1997)
  * [4] Benedikt, M., Watzenig, D., Zehetner, J., Hofer, A.: NEPCE - A nearly energy-preserving coupling element for weak-coupled problems and co-simulation. In: Proceedings of the International Conference on Computational Methods for Coupled Problems in Science and Engineering. pp. 1–12 (Jun 2013)
  * [5] Busch, M.: Performance Improvement of Explicit Co-simulation Methods Through Continuous Extrapolation. In: IUTAM Symposium on solver-coupling and co-simulation. IUTAM Bookseries, vol. 35, pp. 57–80. IUTAM (2019). https://doi.org/10.1007/978-3-030-14883-6_4, iUTAM Symposium on Solver-Coupling and Co-Simulation, Darmstadt, Germany, September 18-20, 2017
  * [6] Éguillon, Y., Lacabanne, B., Tromeur-Dervout, D.: IFOSMONDI: A Generic Co-simulation Approach Combining Iterative Methods for Coupling Constraints and Polynomial Interpolation for Interfaces Smoothness. In: 9th International Conference on Simulation and Modeling Methodologies, Technologies and Applications. pp. 176–186. SCITEPRESS - Science and Technology Publications, Prague, Czech Republic (Jul 2019), https://doi.org/10.5220/0007977701760186
  * [7] Gomes, C., Thule, C., Broman, D., Larsen, P.G., Vangheluwe, H.: Co-simulation: a survey. ACM Computing Surveys (CSUR) 51(3), 1–33 (2018)
  * [8] Gu, B., Asada, H.H.: Co-Simulation of Algebraically Coupled Dynamic Subsystems Without Disclosure of Proprietary Subsystem Models. Journal of Dynamic Systems, Measurement, and Control 126(1), 1–13 (2004). https://doi.org/10.1115/1.1648307
  * [9] Kübler, R., Schiehlen, W.: Two methods of simulator coupling. Mathematical and Computer Modelling of Dynamical Systems 6(2), 93–113 (2000). https://doi.org/10.1076/1387-3954(200006)6:2;1-M;FT093
  * [10] Oosterlee, C.W., Washio, T.: Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows. SIAM Journal on Scientific Computing 21(5), 1670–1690 (2000). https://doi.org/10.1137/S1064827598338093
  * [11] Sadjina, S., Kyllingstad, L.T., Skjong, S., Pedersen, E.: Energy conservation and power bonds in co-simulations: non-iterative adaptive step size control and error estimation. Engineering with Computers 33(3), 607–620 (2017)
  * [12] Schierz, T., Arnold, M., Clauß, C.: Co-simulation with communication step size control in an FMI compatible master algorithm. pp. 205–214 (11 2012). https://doi.org/10.3384/ecp12076205
  * [13] Schweizer, B., Lu, D.: Predictor/corrector co-simulation approaches for solver coupling with algebraic constraints. ZAMM Zeitschrift fur Angewandte Mathematik und Mechanik 95(9), 911–938 (2015). https://doi.org/10.1002/zamm.201300191
  * [14] Sicklinger, S., Belsky, V., Engelman, B., Elmqvist, H., Olsson, H., Wüchner, R., Bletzinger, K.U.: Interface Jacobian-based Co-Simulation. Ph.D. thesis (2014). https://doi.org/10.1002/nme
  * [15] Viel, A.: Implementing stabilized co-simulation of strongly coupled systems using the functional mock-up interface 2.0. In: Proceedings of the 10 th International Modelica Conference. pp. 213–223. Linköping University Electronic Press (March 10 – March 12, Lund, Sweden 2014)

##  Appendix A Parameters of the PETSc non-linear solvers

The JFMs mentionned in this document (see definition in 2.1) refer to PETSc
non-linear solvers, so-called ’SNES’ in the PETSc framework.

The parameters of these methods where the default one, except the explicitely
mentionned ones. The following tables recaps these options. For furthe
definition of their meaning, see [2, 1, 10].

Table 3: Parameters of the NewtonLS method PETSc argument: -snes_linesearch_<...> | Description | Value
---|---|---
`type` |  Select line search type | `bt`
`order` |  Selects the order of the line search for bt | $3$
`norms` |  Turns on/off computation of the norms for basic line search | `TRUE`
`alpha` |  Sets alpha used in determining if reduction in function norm is sufficient | $0.0001$
`maxstep` |  Sets the maximum stepsize the line search will use | $10^{8}$
`minlambda` |  Sets the minimum lambda the line search will tolerate | $10^{-12}$
`damping` |  Damping factor used for basic line search | $1$
`rtol` |  Relative tolerance for iterative line search | $10^{-8}$
`atol` |  Absolute tolerance for iterative line search | $10^{-15}$
`ltol` |  Change in lambda tolerance for iterative line search | $10^{-8}$
`max_it` |  Maximum iterations for iterative line searches | $40$
`keeplambda` |  Use previous lambda as damping | `FALSE`
`precheck_picard` |  Use a correction that sometimes improves convergence of Picard iteration | `FALSE`

Table 4: Parameters of the Anderson method PETSc argument: -snes_anderson_<...> | Description | Value
---|---|---
`m` |  Number of stored previous solutions and residuals | $30$
`beta` |  Anderson mixing parameter | $1$
`restart_type` |  Type of restart | `NONE`
`restart_it` |  Number of iterations of restart conditions before restart | $2$
`restart` |  Number of iterations before periodic restart | $30$
Table 5: Parameters of the Ngmres method (not Ngmres with line search) PETSc argument: -snes_ngmres_<...> | Description | Value
---|---|---
`select_type` |  Choose the select between candidate and combined solution | `DIFFERENCE`
`restart_type` |  Choose the restart conditions | `DIFFERENCE`
`candidate` |  Use NGMRES variant which combines candidate solutions instead of actual solutions | `FALSE`
`approxfunc` |  Linearly approximate the function | `FALSE`
`m` |  Number of stored previous solutions and residuals | $30$
`restart_it` |  Number of iterations the restart conditions hold before restart | $2$
`gammaA` |  Residual tolerance for solution select between the candidate and combination | $2$
`gammaC` |  Residual tolerance for restart | $2$
`epsilonB` |  Difference tolerance between subsequent solutions triggering restart | $0.1$
`deltaB` |  Difference tolerance between residuals triggering restart | $0.9$
`single_reduction` |  Aggregate reductions | `FALSE`
`restart_fm_rise` |  Restart on residual rise from x_M step | `FALSE`
Table 6: Parameters of the Ngmres with linsearch method PETSc argument: -snes_ngmres_<...> | Description | Value
---|---|---
`select_type` |  Choose the select between candidate and combined solution | `LINESEARCH`
⋮
All other options of table A are the same
⋮
PETSc argument: -snes_linesearch_<...> | Description | Value
`type` |  Select line search type | `basic`
`order` |  Selects the order of the line search for bt | $0$
`norms` |  Turns on/off computation of the norms for basic linesearch | `TRUE`
`maxstep` |  Sets the maximum stepsize the line search will use | $10^{8}$
`minlambda` |  Sets the minimum lambda the line search will tolerate | $10^{-12}$
`damping` |  Damping factor used for basic line search | $1$
`rtol` |  Relative tolerance for iterative line search | $10^{-8}$
`atol` |  Absolute tolerance for iterative line search | $10^{-15}$
`ltol` |  Change in lambda tolerance for iterative line search | $10^{-8}$
`max_it` |  Maximum iterations for iterative line searches | $1$
`keeplambda` |  Use previous lambda as damping | `FALSE`
`precheck_picard` |  Use a correction that sometimes improves convergence of Picard iteration | `FALSE`