11institutetext:  Leiden Observatory, University of Leiden, Niels Bohrweg 2,
2333 CA Leiden, the Netherlands
11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>
#  Venice: a multi-scale operator-splitting algorithm for multi-physics
simulations

Maite J. C. Wilhelm 11 S. Portegies Zwart 11

(Received XXXXX; accepted XXXXXX)

###### Abstract

Context. We present Venice, an operator splitting algorithm to integrate a
numerical model on a hierarchy of timescales.

Aims. Venice allows a wide variety of different physical processes operating a
different scales to be coupled on individual and adaptive time-steps. It
therewith mediates the development of complex multi-scale and multi-physics
simulation environments with a wide variety of independent components.

Methods. The coupling between various physical models and scales is dynamic,
and realized through (Strang) operators splitting using adaptive time steps.

Results. We demonstrate the functionality and performance of this algorithm
using astrophysical models of a stellar cluster, first coupling gravitational
dynamics and stellar evolution, then coupling internal gravitational dynamics
with dynamics within a galactic background potential, and finally combining
these models while also introducing dwarf galaxy-like perturbers. These tests
show numerical convergence for decreasing coupling timescales, demonstrate how
Venice can improve the performance of a simulation by shortening coupling
timescales when appropriate, and provide a case study of how Venice can be
used to gradually build up and tune a complex multi-physics model. Although
the examples couple complete numerical models, Venice can also be used to
efficiently solve systems of stiff differential equations.

###### Key Words.:

Methods: numerical — Techniques: miscellaneous

###### doi:

https://doi.org/10.1016/j.procs.2015.05.273

###### doi:

https://doi.org/10.1016/j.jpdc.2012.12.011

###### doi:

10.1088/0067-0049/216/2/29

###### doi:

10.1098/rsta.2013.0378

###### doi:

https://doi.org/10.1016/j.jcp.2015.03.032

###### doi:

https://doi.org/10.1016/j.apm.2016.11.011

###### doi:

10.1051/0004-6361/202038553

###### doi:

10.1093/pasj/59.6.109510.48550/arXiv.0706.2059

###### doi:

10.1051/0004-6361/20142383110.48550/arXiv.1407.7105

###### doi:

https://doi.org/10.1016/j.mcm.2010.09.005

###### doi:

10.1046/j.1365-8711.2001.04022.x

###### doi:

10.1007/BFb0116406

###### doi:

10.1086/147911

###### doi:

10.1086/177173

###### doi:

10.1016/j.newast.2012.05.009

###### doi:

10.1111/j.1365-2966.2011.20137.x

###### doi:

10.1051/0004-6361/201321252

###### doi:

10.48550/arXiv.1808.08028

###### doi:

10.1126/science.aau3206

###### doi:

10.1016/j.newast.2008.10.006

###### doi:

10.1016/j.cpc.2012.09.024

###### doi:

10.1016/j.cnsns.2020.105240

###### doi:

10.1016/j.cnsns.2020.105240

###### doi:

10.2118/182714-MS

###### doi:

10.1137/0705041

###### doi:

10.1093/mnras/224.1.193

###### doi:

10.1051/0004-6361/201218966

###### doi:

10.1051/0004-6361/201834641

###### doi:

10.1103/PhysRev.159.98

###### doi:

10.1177/0037549719881204

##  1 Introduction

Improvements in computational capabilities over recent decades have allowed
numerical models to increase in scale and complexity. This has mostly come in
the form of models covering a range of physical processes – multi-physics
models – and physical scales multi-scale models.

Many multi-physics models use the concept of operator splitting. Instead of
modelling all physical processes in concert, they are solved sequentially.
Operator splitting introduces a numerical error, but this error typically
decreases as the timescale on which the operators are split decreases. If the
operators are applied sequentially for the same time-steps, the error is
linear in the time-steps. Applying the operators for fractional time steps in
interwoven patterns, much like a Verlet scheme Verlet (1967), results in
errors quadratic in the time step (e.g Strang splitting, Strang 1968). Higher-
order operator splitting methods exist, but often in limited cases or non-
trivial forms (Jia and Li 2011; Christlieb et al. 2015; Portegies Zwart et al.
2020a).

Even with traditional Strang splitting, all the operators are applied at
roughly the same timescale, which necessarily must be the shortest common
timescale of the system. In addition, the adopted time step between two
solvers is symmetric. As multi-physics models cover increasingly more physical
processes, they tend to become increasingly more multi-scale as well. In
addition, the range of scales in the problem may dramatically change with
time. In that case, a time step that dynamically varies with time and across
coupling method leads to considerable speed-up and smaller errors.

This issue of a large range of timescales has long been recognised in
modelling systems with a large number of self-gravitating bodies (Aarseth
1985). In response, schemes were developed where particles could be evolved on
their individual, appropriate, time-steps (Wielen 1967; McMillan 1986;
Pelupessy et al. 2012). Additionally, the bridge method (Fujii et al. 2007)
has been developed to couple gravitationally interacting systems of different
physical scales. It splits the combined (Hamiltonian) operator into two
isolated, self-consistent systems and interaction terms. The two systems can
be solved independently, using the preferred methods for each, but they
interact on larger timescales. Although quite powerful, bridge, has several
important shortcomings. These include its fixes and constant coupling time
scale, the rigidity of coupled schemes used, and the fixation of the data set
on which each differential solver operates. As a consequence, the bridge once
constructed, is rigid in the spatial and temporal domain, and in methodology.
This rigidity limits its application in modern simulation endeavors.

With the gradual increase in the sophistication of domain-specific numerical
solvers, it has become increasingly challenging to develop simulation models
incorporating new physics, and employing new algorithms. Expanding those model
also has become more complex, because the higher sophistication renders the
underlying numerical framework more fragile. Without thorough modularization,
we move into a direction in which simulation codes become unwieldingly
complex.

In this paper we introduce Venice, an algorithm to split an arbitrary number
of operators on a wide range of timescales. We explore coupling of operators
on dynamical variable timescales, and discuss how bridge is incorporated into
the framework. We present a reference implementation in the amuse framework,
using it to demonstrate the capabilities and measure the performance of the
algorithm. Venice improves performance in multi-physics simulations by
allowing different components to couple at individually relevant timescales
instead of the shortest common timescale. Possibly even more important than
performance is the possibility to maintain the framework, verify its working
and validate the results. A operational version including research project
using Venice is presented in Huang et al. (2024), where we develope a planet
population synthesis analysis of planet formation in young stellar clusters.
We start, however, with a brief overview of exiting frameworks for addressing
(astro)physical problems with their advantages and disadvantages.

##  2 An brief overview of multi-scale simulation frameworks

Portegies Zwart (2018) argued for a modular approach to multi-scale, multi-
physics simulations. He envisions an environment in which highly optimised
modules address with a limited scope are combined, using a common interface,
into a more complex framework that cover the full palette of physical
processes at their appropriate scale, spatially and temporally. By
incorporating at least two physically motivated solvers per domain, such a
framework would subsequently complies to “Noah’s Ark” philosophy for code
coupling paradigms (Portegies Zwart 2018). Having two or more solvers for
addressing the same physical domain at the same spatial or temporal scale
allows for independent verification and validation of the specific part of the
simulation.

While the whole purpose of performing numerical simulations is to go beyond
what we can test from fundamental principles, an “Noah’s Ark” principle should
be considered one of the fundamental aspects of a simulation framework
environment.

We will briefly discuss a number of such frameworks, which comply in part to
the above criteria. A more extensive overview of such frameworks across
disciplines is presented in Babur et al. (2015). Each of these frameworks has
its advantages and focus area. For example, see Xianmeng et al. (2021) for
simulating power plant reactors, DEM for computational fluid dynamics (Peters
et al. 2018), aluminum electrolysis (Einarsrud et al. 2017), or artery
dynamics in blood flow (Veen and Hoekstra 2020) (discussed below), and the
general implicit coupling framework for multi-physics problems (Rin et al.
2017).

We will briefly discuss the three frameworks that most closely comply to the
listed requirements for a multi-purpose simulation environment. These include
Cactus, the multiscale modeling framework, and amuse.

Cactus (Goodale et al. 2003) is a general purpose problem solving environment
for the physical sciences (Allen et al. 2001). It is modular and supports
parallel computing across the network with a variety of hardware. Cactus had a
last update (version 4.12 of May 2022, see http://www.cactuscode.org). Cactus
offers several advantages over other frameworks. It uses a component-based
architecture, where the core functionality is provided by a central ”flesh”
and various ”thorns” (modules) that can be independently developed, compiled,
and loaded. The thorns can be parallelized, for multi-core architectures and
distributed computing environments. Its high level of abstractions and
interoperability with other tools simplifies the development of complex
simulations,

Closest to the amuse philosophy is probably the Multiscale Modeling and
Simulation Framework (MMSF, Borgdorff et al. 2013; Chopard et al. 2014), which
has been applied in the muscle3 environment (Veen and Hoekstra 2020). However,
in contrast to amuse’s structure of a driver script instructing a set of
models to run and exchange information, it instead manages a set of
distributed, mostly independent models which exchange messages via externally
configurable connections, which block if a result they need is not yet
available.

The Astrophysics Multiphysics Simulation Environment (amuse for short,
Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013) has similar
advantages as Cactus and the MMSF, it has the additional advantage of dynamic
adaptation and a high-level programming interface. Although specifically
designed for applications in astronomy it branches out to oceanography and
hydrology. While amuse provides the capability of coupling an arbitrary number
of codes, it still relies on the coupling of an array of solvers with a
second-order scheme with constant time steps. The amuse framework provides a
first implementation of “Noah’s Ark”, providing a pre-determined communication
protocol with standardized rules and a common data model.

amuse adopts bridge for its code coupling to first, second and any even higher
order (Portegies Zwart et al. 2020a). But it fails on time-dependent
flexibility needed for modern multi-scale and multi-physics simulations. Part
of this problem has been addressed in the nemesis module in amuse (van Elteren
et al. 2019).

Nemesis uses a global structure and sub-structures to couple the stars and
planetary system dynamics through a cascade of Bridge patterns. The global
structure contains stars and planets integrated by a direct $N$-body solver,
while substructures are integrated using a symplectic $N$-body solver coupled
to the global system via a symplectic Bridge Portegies Zwart et al. (2020b).
Nemesis dynamically combines each planetary system’s direct, symplectic and
secular integrators at runtime. Escaping planets are removed from the local
system and incorporated into the global $N$-body integrator, while the
remaining planets continue with their native solver. This is a very dynamic
environment in which codes are spawned and killed at runtime.

With venice we go further into “Noah’s Ark” paradigm, by allowing an operator
based approach to code coupling. By relaxing the coupling strategy, and
allowing for a arbitrary number of codes to be coupled with an arbitrary
number of other codes with arbitrary coupling times makes venice more
flexible, transparant and modular than any of the above mentioned endeavors.

##  3 Algorithm

We here discuss the generalized algorithm for a multi-scale operator approach.

###  3.1 Operator splitting

The Venice algorithm is based on the concept of operator splitting. Consider a
differential equation of the following form:

$\frac{dy}{dt}=L_{1}\left(y\right)+L_{2}\left(y\right).$ (1)

The exact solution to this differential equation is:

$y\left(\Delta t\right)=e^{\left(\hat{L}_{1}+\hat{L}_{2}\right)\Delta
t}y\left(0\right).$ (2)

An approximate solution would be:

$y\left(\Delta t\right)=e^{\hat{L}_{1}\Delta t}e^{\hat{L}_{2}\Delta
t}y\left(0\right)+\mathcal{O}\left(\Delta t\right).$ (3)

This is equivalent to applying the separate operators in sequence. The error
of this scheme is first order in the time-step. Strang (1968) developed a
second order scheme, of the following form:

$y\left(\Delta t\right)=e^{\hat{L}_{1}\Delta t/2}e^{\hat{L}_{2}\Delta
t}e^{\hat{L}_{1}\Delta t/2}y\left(0\right)+\mathcal{O}\left(\Delta
t^{2}\right).$ (4)

In this way one can split two operators. In a multi-physics model, we will
often have more than two terms in an operator. Consider the following
differential equation:

$\frac{dy}{dt}=\sum_{i}l_{i}\left(y\right).$ (5)

We can split this operator by partitioning terms between sub-operators
$\hat{L}_{1}$ and $\hat{L}_{2}$, and recursively splitting sub-operators
consisting of more than one term. This retains the second order scaling
because no terms of $\mathcal{O}\left(\Delta t\right)$ are introduced.

Partitioning between $\hat{L}_{1}$ and $\hat{L}_{2}$ can be done in a variety
of ways, but the most practical is to split off one term into $\hat{L}_{1}$,
and leave the rest in $\hat{L}_{2}$. This minimises the number of operator
applications; if $L_{1}$ ever has to be split again, this has to be done
twice, increasing the number of operators by four. If $L_{2}$ has to be split,
the number of operators only increases by two. For $N$ operators, the split
operator will then be of the following form:

$\begin{split}y\left(\Delta t\right)=&e^{\hat{l}_{1}\Delta
t/2}e^{\hat{l}_{2}\Delta t/2}...e^{\hat{l}_{N-1}\Delta
t/2}e^{\hat{l}_{N}\Delta t}\times\\\ &e^{\hat{l}_{N-1}\Delta
t/2}...e^{\hat{l}_{2}\Delta t/2}e^{\hat{l}_{1}\Delta t/2}.\end{split}$ (6)

###  3.2 Adaptive time operator splitting

Strang splitting also naturally lends itself to adaptive time-stepping
schemes. For this purpose, we separate operators into a slow and a fast set
with respect to a given time-step $\Delta t$; a slow operator can be applied
for $\Delta t$ within some required tolerance, while a fast operator can not.
All slow operators can be partitioned into $\hat{L}_{1}$ and applied for
$\Delta t$, and all fast operators can be partitioned into $\hat{L}_{2}$. They
might not be slow operators at $\Delta t/2$, either, but we can continue
splitting the operators recursively until there are no more fast operators.
Successively halving the time-step results in an exponential decrease, giving
the algorithm a large dynamical range in timescales.

Determining slow and fast operators is not trivial. It is not necessarily
based only on how rapidly an operator works on its section of the data vector,
but also on how other operators work on the same section of the data vector. A
similar problem was considered by Jänes et al. (2014) for Newtonian gravity
n-body models. They pointed out that even if all particles evolve rapidly
(they gave the example of a cluster of hard stellar binaries, which all evolve
on the orbital timescale of the hard binaries), the interactions between stars
in different binaries operate on timescales similar to those of single stars
in the same cluster. They instead proposed an algorithm based on connected
components, and showed that it reduced the wall-clock run-time of a simulation
of a cluster of hard binaries by up to two orders of magnitude, compared to
older adaptive time-stepping schemes.

The algorithm would decompose the particles into subsystems given a set of
particles, a definition of the interaction timescale $\tau_{G,ij}$ between two
particles $i$ and $j$ (e.g. the minimum of the flyby and free-fall times), and
an evolution timescale $\Delta t$. Within a subsystem, or connected component,
each particle $i$ would interact with at least one other particle $j$ with
$\tau_{G,ij}<\Delta t$, while between two subsystems, there were no pairs $ij$
with $\tau_{G,ij}<\Delta t$. This does put particles with $\tau_{G,ij}>\Delta
t$ in the same subsystem, but always because there is some path of particle
pairs with timescales $<\Delta t$ connecting the two.

In the Newtonian gravity n-body model the flyby and free-fall time were
characteristic timescales that could be straightforwardly derived from the
model. For other models this may not be as simple. In Venice, we leave this as
a system parameter. Given $N$ operators, this takes the form of a symmetric
$N\times N$ matrix $\tau$, where $\tau_{ij}=\tau_{ji}$ is the coupling
timescale between operators $\hat{l}_{i}$ and $\hat{l}_{j}$. Two uncoupled
operators (for example, that describe system components that don’t directly
interact), have a coupling timescale of infinity; they can be applied for an
arbitrarily large time-step $\Delta t$ without having to be coupled 111Note
the difference between timescales and time-steps: timescales are a derived
physical quantity of a system, while time-steps are a property of the
algorithm, external to the physics of the system..

Decomposition into connected components results in a set $\hat{C}_{i}$, where
each is a set consisting of $>1$ operator (where operators within a
$\hat{C}_{i}$ interact on timescales $<\Delta t$, but operators in different
$\hat{C}_{i}$’s do not), and a rest set $\hat{R}$ of operators that are not
connected to any other (and do not interact with operators on timescales
$<\Delta t$). Each connected component $\hat{C}_{i}$ is fast compared to
$\Delta t$, and has to be partitioned in $\hat{L}_{1}$. The operators in the
rest set $\hat{R}$ are all slow compared to $\Delta t$, and can be partitioned
into $\hat{L}_{2}$. This results in the following partitioning:

$\frac{dy}{dt}=e^{\left(\sum_{i}\hat{C}_{i}\Delta
t/2\right)}e^{\left(\hat{R}\Delta t\right)}e^{\left(\sum_{i}\hat{C}_{i}\Delta
t/2\right)}.$ (7)

The operators in the rest set $R$, and the connected sets of operators
$C_{1}..C_{K}$, must again be split, which can be done by linear splitting
(Eq. 3) or Strang splitting (Eq. 4). Then, each individual $C_{i}$ can again
be split into connected components. Note that if the connected components
$C_{i}$ are split by Strang splitting, one is split at $\Delta t/2$, while the
others are split at $\Delta t/4$.

In Algorithm 1, we show the basic Venice algorithm. For simplicity we use
first order splitting of the set $C_{1}..C_{K}$ and of $R$.

Algorithm 1 The basic Venice algorithm, with $1^{\textrm{st}}$ order coupling.

1:procedure evolve_cc($L$, $\Delta t$, $\tau$)

2: $\triangleright$ Split codes into connected components

3: $C_{1}$..$C_{K}$, $R$ = split_cc($L$, $\Delta t$, $\tau$)

4: $\triangleright$ Evolve connected components for half time

5: for $C_{i}$ in $C_{1}$..$C_{k}$ do

6: evolve_cc($C_{i}$, $\Delta t/2$, $\tau$)

7: end for

8: $\triangleright$ Evolve rest codes for full time

9: for $l_{i}$ in $R$ do

10: evolve($l_{i}$, $\Delta t$)

11: end for

12: $\triangleright$ Evolve connected components

13: $\triangleright$ for another half time

14: for $C_{i}$ in $C_{1}$..$C_{k}$ do

15: evolve_cc($C_{i}$, $\Delta t/2$, $\tau$)

16: end for

17:end procedure

In Figure 1 we show an illustration of a system of 11 operators, coupled on
different (unit-less for clarity) timescales. This system decomposes into two
subsystems at $\Delta t=2^{-2}\equiv 0.25$, meaning that the two subsystems
exchange information with an interval of $0.25$ (note that time is considered
a dimensionless quantity here). Most of the operators in these subsystems are
then split at a time-step of $\Delta t=2^{-5}\approx 0.03$, except for
operators 6 and 7, which are more tightly coupled and only split at a time-
step of $\Delta t=2^{-8}\approx 0.004$.

Figure 1: Left: a representation of the coupling timescale matrix $\tau$ of an
operator set as an undirected, weighted graph. Edges are labelled with the
coupling timescale; some edges between triangulated nodes share labels for
clarity. The coupling timescale between unconnected nodes is $\infty$. Right:
a tree representation of the recursive splitting of connected components, for
the system on the left. Dotted lines denote skipped layers of nodes with a
single parent and child node.

###  3.3 Kick interactions

The bridge scheme was introduced by Fujii et al. (2007) in order to couple
different gravitational systems, such as stellar clusters within a galaxy. It
uses operator splitting to split the system into individual subsystems and
interaction terms, each with their own time-steps. In this way, a stellar
cluster’s internal dynamics can be solved on a shorter timescale than the
cluster’s orbit in the galaxy, or the dynamics of the galaxy. This method was
expanded by Pelupessy and Portegies Zwart (2012) to also allow gravitational
coupling with smoothed-particle hydrodynamical models. Portegies Zwart et al.
(2020a) expanded the method to higher-order coupling, and have introduced
augmented coupling patterns that allow for external forces (e.g. the YORP
effect and general relativity into an otherwise Newtonian solver). We discuss
how such a splitting can be incorporated into Venice.

A bridged system can be seen as the following operator:

$l_{i}=\sum_{j}\left(l_{i}^{\left(j\right)}+\sum_{k}\tilde{l}_{i}^{\left(jk\right)}\right),$
(8)

where $y_{i}^{\left(j\right)}\cap y_{i}^{\left(k\right)}=0\quad\forall\quad
j\neq k$, i.e., all $l_{i}^{\left(j\right)}$ use and change a different
subvector of $y_{i}$. These are the individual n-body systems. The operators
$\tilde{l}_{i}^{\left(jk\right)}$, however, use both $y_{i}^{\left(j\right)}$
and $y_{i}^{\left(k\right)}$, but affect only $y_{i}^{\left(j\right)}$. These
are the kick interactions that work on n-body system $j$.

If we split the kick interactions $\hat{\tilde{l}}_{i}^{\left(jk\right)}$ into
$L_{1}$ and the evolutions $\hat{l}_{i}^{\left(j\right)}$ into $L_{2}$, we
obtain the following second-order operator splitting (also see Portegies Zwart
et al. 2020a, for higher-order couplings):

$\begin{split}e^{l_{i}\Delta
t}=&\left(\prod_{j}\prod_{k}e^{\hat{\tilde{l}}_{i}^{\left(jk\right)}\Delta
t/2}\right)\left(\prod_{j}e^{\hat{l}_{i}^{\left(j\right)}\Delta
t}\right)\times\\\
&\left(\prod_{j}\prod_{k}e^{\hat{\tilde{l}}_{i}^{\left(jk\right)}\Delta
t/2}\right).\end{split}$ (9)

If an operator equivalent to a bridged system is present, it can be desirable
to split it explicitly for flexibility and transparency. The individual
operators $l_{i}^{\left(j\right)}$ can be treated much in the same way as
other operators $l_{i}$, but the kick operators
$\tilde{l}_{i}^{\left(jk\right)}$ require special treatment. Formally they
must be applied directly before and after the bridged operators, without any
other, non-bridged, operator in between. However, this becomes impossible when
a bridged subsystem is also coupled to another operator on a shorter
timescale222Note the implicit assumption here that the bridge timescale is the
same as the previously defined coupling timescales $\tau_{ij}$. In that case,
the other operator must be included in the second factor of Equation 9, and
that factor must be decomposed in connected components following Equation 7.

The question now is where to apply kick operators. We propose to place them
between the evolution of the connected components and the rest set. This
preserves Strang splitting if two operators are in different connected
components. This is most obvious if both are in the rest set, where two kick
operators flank the evolution operators. If both are in different $C_{i}$’s,
the kick operators are applied without evolution operators in between, which
instead are applied for half a time-step before and after the kick operators.
This preserves the 2nd order bridge scheme.

###  3.4 Dynamic coupling timescales

The decomposition of operators into connected components happens in every
recursive evaluation, including those that happen after some operators have
been (partly) applied. This means that, if the timescale matrix $\tau_{ij}$
were to change as the system evolves, the structure of the connected
components can also change. In this way Venice can create a time-stepping
scheme that adapts to the behaviour of the system. For example, if one
operator represents stellar evolution, this operator can be coupled on long
timescales so long as all stars remain on the main sequence, on shorter
timescales when one or more stars move to more rapidly evolving stages of
evolution, and then again on longer timescales when those stars become inert
stellar remnants.

The coupled system remains consistent (with all operators applied for the same
total time) so long as the connected components within one recursive function
call of evolve_cc (see Algorithm 1) are not changed (i.e. operators $l_{i}$ do
not move between different connected components $C_{i}$ and/or the rest set
$R$). However, between two subsequent recursive function calls of evolve_cc of
a connected component $C_{i}$, the decomposition into further connected
components may be different.

In Algorithm 2, we show the Venice algorithm extended with kick interactions
and dynamic coupling timescales. For clarity, we present the algorithm with
first-order splitting of $R$ and the set $C_{1}..C_{K}$.

Algorithm 2 The extended Venice algorithm, with $1^{\textrm{st}}$ order
coupling.

1:procedure evolve_cc($L$, $\Delta t$, $\tau$)

2: $\triangleright$ Update coupling timescales

3: update_timescales(L, $\Delta t$, $\tau$)

4: $\triangleright$ Split codes into connected components

5: $C_{1}$..$C_{K}$, $R$ = split_cc($L$, $\Delta t$, $\tau$)

6: $\triangleright$ Evolve connected components for half time

7: for $C_{i}$ in $C_{1}$..$C_{k}$ do

8: evolve_cc($C_{i}$, $\Delta t/2$, $\tau$)

9: end for

10: $\triangleright$ Kick between connected components for half time

11: for $l_{i}$ in $l$ do

12: for $l_{j}$ in $l$ do

13: if $l_{i}$ and $l_{j}$ not in same $C_{k}$ then

14: kick($l_{i}$, $l_{j}$, $\Delta t/2$)

15: end if

16: end for

17: end for

18: $\triangleright$ Evolve rest codes for full time

19: for $l_{i}$ in $R$ do

20: evolve($l_{i}$, $\Delta t$)

21: end for

22: $\triangleright$ Kick between connected components for another half time

23: for $l_{i}$ in $l$ do

24: for $l_{j}$ in $l$ do

25: if $l_{i}$ and $l_{j}$ not in same $C_{k}$ then

26: kick($l_{i}$, $l_{j}$, $\Delta t/2$)

27: end if

28: end for

29: end for

30: $\triangleright$ Evolve connected components for another half time

31: for $C_{i}$ in $C_{1}$..$C_{k}$ do

32: evolve_cc($C_{i}$, $\Delta t/2$, $\tau$)

33: end for

34:end procedure

###  3.5 amuse implementation

Thus far we have discussed Venice as an operator splitting algorithm. However,
it can also be used as a coupling scheme for physical models that may not
strictly solve a proper partial differential equation. The amuse platform
provides a standardised structure for (interfaces to) physical models, making
it an excellent environment to develop a reference implementation of the
algorithm.

In the previous sections, we tacitly assumed that all operators affect the
same data vector. However, in amuse every code has its own data vector
(typically in the form of particle or grid data). Coupling must be done
explicitly, by copying data from one code to another. In this implementation
we use the channels in amuse. These are defined from one data structure to
another, and can copy any subset of data from one structure to the other. They
must be explicitly called, which we do between operations. We could have
chosen a lazy approach, where before the evolution of a code we copy all
channels to that code. However, this results in an end state where the data of
different codes is different, and it is non-trivial to reconstruct which code
contains the most recent state. Instead, we take an active approach, where we
copy all channels from a code after its evolution step. This can introduce
unnecessary copy operations, but is less prone to mistakes.

In many situations the structure of a code’s data set will be updated
throughout the simulation (e.g. adaptive mesh refinement codes, particle-based
stellar winds and sink particles). To resolve this, every pair of codes can
have a function that synchronises the data sets, defined by the user. This is
done before channels are copied to ensure the data structures are compatible.

We also implemented a data output framework. This allows for analysis of fast
subsystems without needing the same short time-step for the entire system (or
the construction of a wrapper class whose only purpose besides evolving is
writing data). We write three types of output: checkpoint files, plot files,
and debug files. Checkpoint files are saved after the evolution of all codes,
and contain the state of the entire system at a synchronised time. Plot files
are saved after the evolution of an entire connected component. To prevent
redundant files, saving checkpoint and plot files is skipped if there is only
a single connected component, consisting of more than one code; in that case,
the state at the end of the function is the same as at the end of the second
recursive call. Debug files are saved after the evolution of individual codes.

The amuse implementation is available on Github:
https://github.com/amusecode/venice/.

##  4 Testing Venice

In this section we present a number of tests of Venice using the amuse
implementation. The scripts to run these tests, and produce the figures, are
also available on the Venice Github page.

###  4.1 Gravity and stellar evolution

Coupled gravity-stellar evolution simulations are one of the oldest types of
astrophysical multi-physics simulations, going back to Terlevich (1987). The
coupling of these systems is relatively straightforward. Assuming a
collisionless system and no changes in momentum during stellar evolution (e.g.
supernova kicks), the only coupling is through the mass. Stars lose mass in
stellar winds and at the end of their lifetime, and this changes the
gravitational acceleration they induce in other stars. Newton’s equations of
motion can account for changing mass in differential form, but most
astrophysical gravitational dynamics models do not include this. Coupling is
typically done through operator splitting.

We use the ph4 (McMillan et al. 2012) gravitational dynamics code and the SeBa
parameterised stellar evolution code (Portegies Zwart and Verbunt 1996; Toonen
et al. 2012). Coupling is done through a channel copying mass from the stellar
evolution code to the gravitational dynamics code.

In Section 4.1.1 we use a constant coupling timescale. This becomes
problematic if the stellar evolution timescale changes. This happens, for
example, when a star leaves the main sequence. In Section 4.1.2, we expand the
model with a dynamic coupling timescale.

For the initial conditions we adopt a 100 star Plummer sphere in virial
equilibrium. Stellar masses are distributed following a Kroupa IMF (Kroupa
2001), between 0.08 and 8 M⊙. We manually set a randomly selected star to a
higher mass, 16 M⊙ in Section 4.1.1 and 50 M⊙ in Section 4.1.2. The Plummer
radius is set such that the free-fall time is 6 Myr (resulting in a Plummer
radius of $\sim$1.25 pc with the 16 M⊙ star, and $\sim$1.47 pc with the 50 M⊙
star). This is also the total run-time of the tests. During this time, the 16
M⊙ star remains on the main sequence, and loses relatively little mass,
allowing us to use relatively long coupling timescales compared to the total
simulation time. The 50 M⊙ star leaves the main sequence and undergoes a
supernova explosion during this time frame. This results in more rapid changes
in mass, requiring very short coupling timescales at certain epochs. In
Section 4.1.2, we demonstrate a timescale update function that can shorten the
coupling timescale appropriately.

####  4.1.1 Constant coupling timescale

In Algorithm 3 we show how a Venice system of gravitational dynamics and
stellar evolution coupled at a constant timescale is set up.

Algorithm 3 Venice setup of gravitational dynamics and stellar evolution
coupled at a constant timescale

1:procedure setup_gravity_stellar_constant( gravity, stellar, dt)

2: system = Venice()

3: system.add_code(gravity)

4: system.add_code(stellar)

5: system.timescale_matrix[0][1] = dt

6: $\triangleright$ Copy mass from code 1 (stellar)

7: $\triangleright$ to 0 (gravity)

8: system.add_channel(1, 0, [”mass”])

9: return system

10:end procedure

In Figure 2 we show three views of the median position error, coupling
timescale, and wall-clock run-time as functions of each other, for both 1st
and 2nd order splitting of the operator rest set $R$ (the splitting order of
the connected components $C_{1}..C_{K}$ does not affect the splitting
procedure for $<4$ operators). We also show the time evolution of the median
position error. The position error is defined as the distance between
initially identical particles in runs with different timescales. We compute
position errors between subsequent values of the coupling timescale, indicated
by the end points of error bars.

Most gravitational systems of more than two bodies are chaotic (Miller 1964).
As a result, any perturbation we introduce in the system triggers an
exponential growth of that perturbation. Changing the timescale from one run
to another is such a perturbation, and therefore triggers exponential growth.
We observe that the median position error grows roughly linearly with time
(although this may be the first-order behaviour of exponential growth at times
much smaller than the Lyapunov time), and that there is a systematic offset
between curves associated with subsequent coupling timescales. Therefore, we
argue that the systematic change is caused by the coupling timescale, and not
by the chaotic behaviour of the system.

The median position error scales approximately linearly with the coupling
timescale, for both 1st and 2nd order coupling. However, at the same time-step
the error of the $2^{\textrm{nd}}$ order method is an order of magnitude
smaller. This comes at the cost of slightly longer wall-clock run-time333The
increase in run-time is due to the gravity model running twice as long. but
the median position error at a given run-time is typically lower for
$2^{\textrm{nd}}$ order than for $1^{\textrm{st}}$ order. In this model
configuration (where the gravity model is the first element of the operator
list $L$), it is the gravity model that is evolved in two steps due to Strang
splitting. This can be explained by the gravity model’s preferred internal
time-step being larger than the largest coupling timescale, so that it
performs a single numerical step every time Venice evolves the model. In this
case, the gravity model does not dominate computation time, so the increase in
run-time is small. This implies that $2^{\textrm{nd}}$ order coupling results
in the best trade-off between run-time and performance unless both models have
comparable run-time and have to shorten their internal time-step.

Figure 2: Convergence and performance of a gravitational dynamics and a
stellar evolution model coupled at a constant timescale. The top panels and
bottom left panel show three views of the median position error, coupling
timescale, and wall-clock run-time as functions of each other. The bottom
right panel shows the time evolution of the median position error in all runs.
The black dashed lines in the top left plot show linear and quadratic scalings
of the error with the timescale.

The order in which operators are split affects the outcome of numerical
integration. To quantify this effect we compare the model above with a model
that is identical, except that the stellar evolution model is the first
operator in $L$. This is done by reversing the order of add_code calls in
Algorithm 3, and results in the stellar evolution model being evolved twice
per Strang splitting, instead of the gravity model. In Figure 3, we show how
the median position error and the wall-clock run-time depend on the coupling
timescale, for the two reversed models, at the same coupling timescale, for
both 1st and 2nd order coupling.

We observe globally the same behaviour as for the convergence with timescale;
the median position errors converges linearly with the timescale, for both 1st
and 2nd order coupling, although 2nd order coupling results in smaller errors.
The magnitude of the errors are also roughly the same at a given coupling
timescale. Finally, we again see that the code being evolved twice has a
longer run-time444If the system is evolved for multiple time-steps $\Delta t$,
there may be two sequential half-time-steps, e.g.: $e^{\hat{L}_{1}\Delta
t/2}e^{\hat{L}_{2}\Delta t}\left(e^{\hat{L}_{1}\Delta t/2}e^{\hat{L}_{1}\Delta
t/2}\right)e^{\hat{L}_{2}\Delta t}e^{\hat{L}_{1}\Delta t/2}$. These could be
combined in a single evolution step, with a potential for faster evolution,
but for larger systems there could be operators in between, and the coupling
timescales may change in between. Therefore such recombination is not
possible.. This shows that the order of operators impacts the efficiency of
the model, and while it also impacts the results, this error converges with
shorter coupling timescales.

Figure 3: Convergence and performance of a gravitational dynamics and a
stellar evolution model coupled at a constant timescale, with reversed
operators. We show, as a function of coupling timescale, the median position
error (top) and wall-clock run-time (bottom). The black dashed lines in the
top plot show linear and quadratic scalings of the error with the timescale.

####  4.1.2 Adaptive coupling timescale

In the test above, we were able to achieve a numerically converging solution
at relatively long coupling timescales because no star evolved off the main
sequence during the integration time. Off the main sequence, stars evolve
rapidly, ultimately going supernova on a stellar free-fall time, or shedding
their envelopes as planetary nebulae. Evolving the entire system at the
shortest possible timescale at which the stellar population can evolve would
result in a large overhead, leading to many unnecessary couplings. Venice can
easily resolve this using an update_timescale function. Setting this to some
fraction $\eta$ of the evolution timescale of the fastest evolving star allows
the system to couple more frequently when needed.

In principle, the most efficient integration would be achieved by coupling
every star on its own timescale; the scheme we described evolves all stars,
including those on the main sequence, on the timescale of the star with the
shortest time-step. However, this increases the overhead of Venice (note that
the connected component algorithm scales quadratically with the number of
operators). The optimal solution is to separate stars into a number of mass
bins, but calibrating the optimal method of distribution is beyond the scope
of this work.

In Algorithm 4 we show how a Venice system of gravitational dynamics and
stellar evolution coupled at a dynamic timescale is set up.

Algorithm 4 Venice setup of gravitational dynamics and stellar evolution
coupled at an adaptive timescale

1:procedure setup_gravity_stellar_adaptive( gravity, stellar, $\eta$)

2: system = Venice()

3: system.add_code(gravity)

4: system.add_code(stellar)

5: $\triangleright$ lambda defines a function, as in Python

6: system.update_timescale[0][1] =

7: lambda code1, code2, dt:

8: $\eta$ * code2.stellar_timesteps.min()

9: $\triangleright$ Copy mass from code 1 (stellar)

10: $\triangleright$ to 0 (gravity)

11: system.add_channel(1, 0, [”mass”])

12: return system

13:end procedure

We use the same initial conditions as in the previous section, except with a
50 M⊙ star instead of a 16 M⊙ star, and rescaling the cluster radius and
stellar velocities such that the dynamical time is still 6 Myr.

In Figure 4 we show the median position error (top) and the wall-clock run-
time (bottom) as a function of the coupling timescale parameter $\eta$, for
both 1st and 2nd order coupling. The median position error scales roughly
linearly for both 1st and 2nd order coupling, but the error does not
monotonically scale with the timescale parameter. The increase in run-time at
$2^{\mathrm{nd}}$ order coupling is again due to gravity being split into two
evolution steps.

Figure 4: Convergence and performance of a gravitational dynamics and a
stellar evolution model coupled at an adaptive timescale. We show, as a
function of the timescale parameter $\eta$, the median position error (top)
and wall-clock run-time (bottom). The median position error is computed for
two subsequent values of the coupling timescale, indicated by the endpoints of
the error bars. The black dashed lines in the top plot show linear and
quadratic scalings of the error with the timescale.

In Figure 5 we show how the coupling timescale evolves as a function of time
for a number of values of $\eta$ (top), and the stellar type of the 50 M⊙ star
(bottom) over the same time interval. Over the first 4 Myr, the star is on the
main sequence and the coupling timescale is mostly constant. At a number of
transitions between phases (before and after the Hertzsprung gap and giant
branch naked helium star phases, sometimes called “stripped stars”, these
days), the coupling timescale drops by up to seven orders of magnitude, but
also quickly increase by orders of magnitude afterwards. Then, when the star
is a black hole, the coupling timescale becomes the same for values of
$\eta\geq 0.1$. This indicates that the coupling timescale dictated by the
stellar evolution model is large compared to the remaining integration time.

Figure 5: Top: The dynamic coupling timescale between the gravitational
dynamics and stellar evolution models as a function of time, for a number of
timescale parameters $\eta$. Bottom: The stellar type through time of a 50 M⊙
star, which is the most massive one in the system of the top panel.

We note that for $\eta>1$, there is a tendency for the model to ‘overshoot’
phases of stellar evolution. Ideally, the ratio of coupling timescales for two
values of $\eta$ should always be about $\eta$ itself, but during e.g. the
transition from the Hertzsprung Gap to Core Helium burning for $\eta=10$, and
the transition from Giant Branch Naked Helium Star to Black Hole for
$\eta=3,10$, the sharp drop to coupling timescales $\ll$1 kyr does not occur,
while for higher values of $\eta$ it does. In these cases, the stellar
evolution model evolves through the rapid evolution phase without recomputing
the coupling timescale. This implies that values of $\eta\lesssim 1$ should be
used.

###  4.2 Bridged gravity

We demonstrate the bridge algorithm in Venice by simulating the dynamics of a
star cluster in a galaxy.

The cluster consists of 100 stars, with masses drawn from a Kroupa IMF between
0.08 and 8 M⊙, distributed following a Plummer distribution with a radius of
10 pc. We again use ph4 to model the cluster’s gravitational dynamics.

We couple this cluster to the analytic Bovy (2015) Milky-Way-like galaxy
potential, which consists of a bulge (power-law, exponentially cut-off density
profile), disk (Miyamoto-Nagai density profile, Miyamoto and Nagai 1975), and
halo (NFW density profile, Navarro et al. 1996).

We initialise the cluster at 8 kpc from the centre of the galaxy, with a
velocity perpendicular to the direction to the centre. We parameterise this
velocity in terms of an effective eccentricity: the velocity is that of a
Keplerian orbit, at apocenter, for a central mass equal to the total mass
within a circular orbit, given the effective eccentricity. If this
eccentricity is 0, the orbit will be circular, but larger effective
eccentricities do not result in elliptic orbits (because the mass enclosed in
the orbit is not constant). We evolve the system for $223$ Myr, one period of
the adopted orbit.

####  4.2.1 Classical bridge

In Algorithm 5 we show how a Venice system of two bridged gravitational
dynamics models, coupled at a constant timescale, is set up. In this case,
code 0 is the galaxy potential, and code 1 is the cluster gravity.

We set the cluster’s orbital velocity such that the effective eccentricity is
0. In these tests, we use the cluster’s centre of mass to investigate the
system’s convergence.

Algorithm 5 Venice setup of two gravitational systems, where code 0 kicks code
1 through bridge, at a constant timescale.

1:procedure setup_bridged_gravity_constant( gravity1, gravity2, dt)

2: system = Venice()

3: system.add_code(gravity1)

4: system.add_code(gravity2)

5: system.timescale_matrix[0][1] = dt

6: $\triangleright$ Kick interactions from code 0 to code 1

7: system.kick[0][1] = dynamic_kick

8: return system

9:end procedure

In Figure 6 we show how the centre of mass error (top) and wall-clock run-time
(bottom) depend on the coupling timescale. Because the only coupling is done
through kick interactions (and there is no shared data between the gravity
codes), the splitting order of the operator rest set $R$ should not affect the
result. For this reason, we only show the result of $1^{\textrm{st}}$ order
splitting. However, we note that by construction the bridge scheme has
$2^{\textrm{nd}}$ order error behaviour.

The position error indeed scales quadratically with the coupling timescale.
Towards long coupling timescales, the run-time of the cluster gravity model
appears to move towards an asymptote. This is the point at which the code no
longer requires shortened time-steps to fit within the coupling timescale, but
can take large internal time-steps.

Figure 6: Centre of mass error (top) and wall-clock run-time (bottom) as a
function of constant coupling timescale, for a star cluster in a background
galaxy potential. Centre of mass errors are computed between subsequent values
of the coupling timescale, indicated by the ends of the error bars. The dashed
black line indicates quadratic scaling of the centre of mass error with the
coupling timescale.

####  4.2.2 Adaptive bridge

Analogous to the time-step in the general n-body problem, the bridge time-step
can potentially cover a large range. If a star cluster’s orbit within a galaxy
is not circular, the appropriate bridge time-step can be smaller close to the
galactic centre than farther away.

Algorithm 6 Venice setup of two gravitational systems, where system 1 kicks
system 2 through bridge, at an adaptive timescale.

1:procedure setup_bridged_gravity_adaptive( gravity1, gravity2, $\eta$)

2: system = Venice()

3: system.add_code(gravity1)

4: system.add_code(gravity2)

5: system.update_timescale[0][1] =

6: lambda code1, code2, dt:

7: $\eta$ * min_freefall_time(code1, code2)

8: $\triangleright$ Kick interactions from code 0 to code 1

9: system.kick[0][1] = dynamic_kick

10: return system

11:end procedure

We base the adaptive coupling timescale on the freefall times of the particles
in one model in the potential of the other model. The freefall time between
two particles $i$ and $j$ is $\sqrt{r_{ij}/a_{ij}}$, where $r_{ij}$ is the
distance between the particles, and $a_{ij}$ is the gravitational
acceleration. The acceleration is computed for the bridge interaction, and the
distance is that to the centre of the central potential. Similar to the
dynamic stellar coupling timescale, we now take the minimum among all
particles and apply a scaling factor $\eta$. We use the initial conditions
described above, and set the cluster’s velocity such that the effective
eccentricity is 0.5.

In Figure 7 we show how the centre of mass error (top) and wall-clock run-time
(bottom) depend on the coupling timescale parameter $\eta$. Similarly to the
case of the constant coupling timescale, the centre of mass error converges
quadratically with the coupling timescale parameter, and the cluster gravity
run-time no longer decreases at large coupling timescales. We note that the
centre of mass error between $\eta=0.3$ and $\eta=1$ is larger than the
initial distance from the galactic centre, while between $\eta=0.1$ and
$\eta=0.3$ it is almost two orders of magnitude smaller than the initial
distance. This implies that a value $\eta<1$ should be used.

Figure 7: Centre of mass error (top) and wall-clock run-time (bottom) as a
function of timescale parameter $\eta$, for a star cluster in a background
galaxy potential. Centre of mass errors are computed between subsequent values
of the coupling timescale, indicated by the ends of the error bars. The dashed
black line indicates quadratic scaling of the centre of mass error with the
coupling timescale.

###  4.3 Combined model

Having demonstrated Venice’s capabilities in small, controlled experiments, we
now demonstrate its ability to build up complex simulations. In this section
we couple the models described above. We aim to simulate the dynamics of a
star cluster in the potential of a galaxy, while following stellar mass loss
due to stellar evolution. We set up the interactions between the cluster’s
gravitational dynamics and stellar evolution models adaptively, in the same
way as in Section 4.1.2, and the interaction between the cluster’s
gravitational dynamics and the background potential as in Section 4.2.2. We
use a stellar coupling timescale parameter ($\eta_{s}$) of 0.3, and a
gravitational coupling timescale parameter ($\eta_{g}$) of 0.03.

In addition to combining the two models, we add two massive particles just
outside the galactic disk, comparable in mass and location to the Large and
Small Magellanic Clouds. We add these in their own gravity code, which is also
coupled to the galaxy, and also perturbs the cluster (the cluster has no
effect on other components). This separation into two codes allows us to
couple the galaxy and these perturbers to the cluster on different timescales.
In Algorithm 7, we show how this systems is set up.

For the cluster, we use the same 100 star, 10 pc Plummer sphere, with masses
from a Kroupa IMF up to 8 M⊙, except that we give three stars masses of 16,
30, and 50 M⊙. The cluster is initially on a circular orbit at 5 kpc from the
galactic centre. The two perturbers have a mass of $10^{11}$ M⊙ each, and are
initially at 50 and 60 kpc circular orbits, 180∘ and 90∘ ahead of the cluster,
respectively. We evolve the system for 5 Gyr.

Algorithm 7 Venice setup of a cluster with stellar evolution in a galaxy
perturbed by dwarf galaxies.

1:procedure setup_bridged_gravity(gravity_cluster, stellar, gravity_galaxy,
gravity_perturbers, $\eta_{s}$, $\eta_{g}$)

2: system = Venice()

3: system.add_code(gravity_cluster)

4: system.add_code(stellar)

5: system.add_code(gravity_galaxy)

6: system.add_code(gravity_perturbers)

7: system.add_channel(1, 0, [”mass”])

8: $\triangleright$ All gravity systems are coupled; stellar evolution only to
cluster gravity

9: system.update_timescale[0][1] = lambda code1, code2, dt:
$\eta_{s}$*code2.stellar_timesteps.min()

10: system.update_timescale[0][2] = lambda code1, code2, dt:
$\eta_{g}$*min_freefall_time(code1, code2)

11: system.update_timescale[0][3] = lambda code1, code2, dt:
$\eta_{g}$*min_freefall_time(code1, code2)

12: system.update_timescale[2][3] = lambda code1, code2, dt:
$\eta_{g}$*min_freefall_time(code1, code2)

13: $\triangleright$ Galaxy and perturbers both kick cluster

14: system.kick[2][0] = dynamic_kick

15: system.kick[3][0] = dynamic_kick

16: $\triangleright$ Galaxy kicks perturbers

17: system.kick[2][3] = dynamic_kick

18: return system

19:end procedure

In Figure 8, we show projections of this system on the galactic plane. The
perturbers interact before having completed a full orbit, kicking the outer
body inwards (towards the inner 10-15 kpc, still outside the cluster orbit),
and kicking the inner body outwards. However, the cluster is still clearly
perturbed; in runs without perturbers, the cluster disperses along a segment
of a circle still 5 kpc in radius, while in runs with perturbers, stars are
scattered throughout the galaxy. Overdensities remain in the inner 5 kpc, and
around the final position of the perturber that was scattered inwards. A
number of former cluster members appear to have become bound to this
perturber.

Figure 8: Maps of the combined simulation projected on the galactic plane. The
left panel shows a wide view, while the right panel is zoomed in on the
centre. The black dashed and dotted lines show the trajectories of the two
perturbers. Points show the positions at 5 Gyr of the cluster stars, where the
colour indicates the model setup (red and grey are with perturbers, orange and
grey are with stellar evolution).

In Figure 9, we show the evolution of the star cluster by following the
standard deviation of the cluster members’ polar coordinates through time.

Figure 9: The standard deviation of the distance to the galactic centre
normalised to the mean distance (top) and of the azimuthal angle in the
galactic plane (bottom) of the cluster members as a function of time. Colour
indicates the model setup (red and grey are with perturbers, orange and grey
are with stellar evolution).

If there are no perturbers, the radial standard deviation oscillates but
remains below 1%, meaning that the stars remain close to their initial 5 kpc
circular orbit. The azimuthal standard deviation increases mostly linearly,
meaning that the cluster steadily disperses along its orbit. After 4 Gyr there
are steps in this standard deviation. This is likely a result from phase
wrapping; if a considerable fraction of stars has angular coordinates around
$-\pi$ or $\pi$, the standard deviation may be biased. Some sections continue
the initial linear trend, so the cluster continues to disperse. Both the
radial and azimuthal standard deviations are smaller when stellar evolution is
included. This is a counterintuitive result, because stellar evolution results
in an increasingly unbound cluster, which should promote dispersal into the
galaxy.

When perturbers are added, the radial and azimuthal standard deviations
increase and become less regular. After 3 Gyr, the radial standard deviation
becomes of order unity, implying that the cluster is entirely scattered
throughout the galaxy. The azimuthal standard deviation begins to deviate from
the case without perturbers at about 2.5 Gyr. Note that the standard deviation
of a uniform distribution between $\left[0,2\pi\right]$ is
$\frac{\pi}{\sqrt{3}}\approx 1.81$, comparable to the value after 2.5 Gyr.

In Figure 10, we show the evolution time-steps of cluster gravity, perturber
gravity, and stellar evolution codes through time, for the run with both
perturbers and stellar evolution. The stellar time-steps drop below 1 yr
during phases of rapid stellar evolution (as seen before in Fig. 5), but can
also be as large as $\sim$100 Myr. The perturber and cluster time-steps are
mostly constant before 3 Gyr (except when the stellar time-steps are small),
reflecting that they mostly interact with the galaxy rather than each other.
After 3 Gyr, these time-steps become much shorter, and this coincides with the
moment the cluster’s radial standard deviation reaches unity. This confirms
that the cluster’s disruption was due to interactions with a perturber.

Figure 10: The evolution timescales of different components of the combined
model (including both stellar evolution and perturbers) through time. The
dashed line indicates the cluster gravity, the dash-dotted line the perturber
gravity, and the dotted line stellar evolution.

##  5 Discussion

###  5.1 Future applications of the amuse implementation

Although almost any amuse simulation could be refactored using Venice, this
does not necessarily yield a more efficient code or a clearer script. Venice
is especially useful in cases with many models or a wide (and dynamic) range
of timescales. The adaptive coupling between gravity and stellar evolution
could be worked out further, for example by finding the ideal distribution of
stars among multiple stellar evolution codes depending on their lifetime
distribution. This method could also be used for simulating the dynamics of
smaller systems including planetary systems, such as the consequences of a
planetary system once the parent star turns into a compact object, or the
evolution of massive stellar triples (and higher-order multiples) in stellar
clusters.

Venice’s modular expandable nature can be used to gradually and hierarchically
build up complex models, allowing both modules and interactions between
modules to be characterised individually. A simple planet formation model
where a single planet migrates and accretes from a static (gas, pebble, or
planetesimal) disk can slowly be expanded to resemble e.g. the Bern model (see
Figure 2 in Emsenhuber et al. 2021, which has 14 physical modules), where the
effect of each module can be characterised individually. Codes can also be
exchanged for more efficient or advanced alternatives, which should not affect
interactions between other parts of the model that do not interact with that
code.

###  5.2 Steps towards a differential equation solver

We introduce Venice; a coupling paradigm between numerical models. The
designed method is also suitable for being applied to differential equation
solvers. The main components that have not been developed yet are the
integration of the individual operators (corresponding with the evolve_model
functions) and the calculation of the coupling timescales $\tau_{ij}$. Here we
sketch how such a method could be constructed.

The individual integration method could be practically any method with time-
step control (e.g. Euler (forward or backward), (embedded) Runge-Kutta, or
Bulirsch-Stoer), although a single method must be chosen for all operators.
The coupling timescale of an operator pair will be the next time-step
(predicted by the chosen method) for the differential equation where the
right-hand side consists of just those two operators. If these operators are
not coupled, this predicted timescale will be determined by the slowest of the
two operators. This does not reflect their actual coupling timescale, which is
still infinite. This was not an issue in the amuse implementation, where the
internal evolution time steps were decoupled from the coupling timescale.

The Venice algorithm does not eliminate stiff parts of a set of differential
equations, but decouples them from non-stiff parts. The integration method
should therefore still be implicit in order to maintain stability. A
(Bulirsch-Stoer-like) extrapolation method may be useful, because it may allow
the slower operator of a fast pair to converge with few iterations, whereas an
Euler/Runge-Kutta method will force the slow operator to take the same short
time-steps as the fast operator.

A particular application for a Venice-based differential equation solver would
be a time-dependent chemical network integrator. This consists of many
operators (chemical reactions) that operate on a range of timescales depending
on reaction rates and chemical abundances.

##  6 Conclusion

We introduced Venice, an operator splitting algorithm for multi-scale, multi-
physics numerical problems, and demonstrated its working with examples from
astrophysics. The simple examples show convergence with decreasing coupling
times, although even with $2^{\mathrm{nd}}$ order (Strang) coupling this
convergence is potentially only $1^{\mathrm{st}}$ order, although
$2^{\mathrm{nd}}$ order coupling typically results in smaller errors (even at
a given wall-clock run-time). We have also shown that updating coupling
timescales can result in appropriate adaptive time-stepping, and that
gravitational bridge interactions can be included as a coupling method. We
then demonstrated how a complex model could be built up with Venice.

The amuse implementation of Venice can now be used as a framework for
simulations coupling a large number of numerical models on a wide range of
timescales. We have also illustrated how Venice can be used to construct a
differential equation solver. This requires additional development, such as
how to treat the timescales of individual evolution and coupling, and which
integration method to use.

## Acknowledgments

We thank Lourens Veen for the useful discussion on the differences and
similarities between muscle3 and venice, and Shuo Huang for early user
feedback. MW is supported by NOVA under project number 10.2.5.12.

The simulations in this work have been carried out on a work station equipped
with an Intel Xeon W-2133 CPU, which consumes 140 W. At most, one out of six
cores was used per simulation, putting power usage at 23 W. The amount of
computer time used to develop and test Venice and to produce the figures in
this work were not recorded, but will be on the order of $\sim$300 CPU hours.
This puts the total power consumption at $\sim$6 kWh. At 0.649 kWh/kg of CO2
emitted (Dutch norm for grey electricity), this comes to $\sim$10 kg of CO2,
comparable to a week of daily commutes.

## References

  * Aarseth (1985) Aarseth, S.J., 1985.  Direct methods for N-body simulations., in: Multiple time scales, p. 377 - 418, pp. 377–418.
  * Allen et al. (2001) Allen, G., Benger, W., Dramlitsch, T., Goodale, T., Hege, H.C., Lanfermann, G., Merzky, A., Radke, T., Seidel, E., 2001.  Cactus grid computing: Review of current development, in: Euro-Par 2001 Parallel Processing: 7th International Euro-Par Conference Manchester, UK, August 28–31, 2001 Proceedings 7, Springer. pp. 817–824.
  * Babur et al. (2015) Babur, O., Smilauer, V., Verhoeff, T., van den Brand, M., 2015\.  A survey of open source multiphysics frameworks in engineering, in: Procedia Computer ScienceVolume 51Issue C, Technische Universiteit Eindhoven. p. 1088–1097.  doi:.
  * Borgdorff et al. (2013) Borgdorff, J., Falcone, J.L., Lorenz, E., Bona-Casas, C., Chopard, B., Hoekstra, A.G., 2013\.  Foundations of distributed multiscale computing: Formalization, specification, and analysis.  Journal of Parallel and Distributed Computing 73, 465–483.  URL: https://www.sciencedirect.com/science/article/pii/S0743731512002924, doi:.
  * Bovy (2015) Bovy, J., 2015.  galpy: A python Library for Galactic Dynamics.  ApJS 216, 29.  doi:, arXiv:1412.3451.
  * Chopard et al. (2014) Chopard, B., Borgdorff, J., Hoekstra, A., 2014.  A framework for multi-scale modelling.  Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 372.  doi:.
  * Christlieb et al. (2015) Christlieb, A.J., Liu, Y., Xu, Z., 2015.  High order operator splitting methods based on an integral deferred correction framework.  Journal of Computational Physics 294, 224–242.  URL: https://www.sciencedirect.com/science/article/pii/S0021999115001795, doi:.
  * Einarsrud et al. (2017) Einarsrud, K.E., Eick, I., Bai, W., Feng, Y., Hua, J., Witt, P.J., 2017.  Towards a coupled multi-scale, multi-physics simulation framework for aluminium electrolysis.  Applied Mathematical Modelling 44, 3–24.  URL: https://www.sciencedirect.com/science/article/pii/S0307904X16306060, doi:.
  * Emsenhuber et al. (2021) Emsenhuber, A., Mordasini, C., Burn, R., Alibert, Y., Benz, W., Asphaug, E., 2021\.  The New Generation Planetary Population Synthesis (NGPPS). I. Bern global model of planet formation and evolution, model tests, and emerging planetary systems.  A&A 656, A69.  doi:, arXiv:2007.05561.
  * Fujii et al. (2007) Fujii, M., Iwasawa, M., Funato, Y., Makino, J., 2007\.  BRIDGE: A Direct-Tree Hybrid N-Body Algorithm for Fully Self-Consistent Simulations of Star Clusters and Their Parent Galaxies.  Publ. Astr. Soc. Japan 59, 1095.  doi:, arXiv:0706.2059.
  * Goodale et al. (2003) Goodale, T., Allen, G., Lanfermann, G., Massó, J., Radke, T., Seidel, E., Shalf, J., 2003.  The cactus framework and toolkit: Design and applications, in: Palma, J.M.L.M., Sousa, A.A., Dongarra, J., Hernández, V. (Eds.), High Performance Computing for Computational Science — VECPAR 2002, Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 197–227.
  * Huang et al. (2024) Huang, S., Portegies Zwart, S., Wilhelm, M., 2024.  The formation and evolution of planetary systems in dense stallar clusters.  submitted to A&A .
  * Jänes et al. (2014) Jänes, J., Pelupessy, I., Portegies Zwart, S., 2014.  A connected component-based method for efficiently integrating multi-scale N-body systems.  A&A 570, A20.  doi:, arXiv:1407.7105.
  * Jia and Li (2011) Jia, H., Li, K., 2011.  A third accurate operator splitting method.  Mathematical and Computer Modelling 53, 387–396.  URL: https://www.sciencedirect.com/science/article/pii/S089571771000436X, doi:.
  * Kroupa (2001) Kroupa, P., 2001.  On the variation of the initial mass function.  MNRAS 322, 231–246.  doi:, arXiv:astro-ph/0009005.
  * McMillan et al. (2012) McMillan, S., Portegies Zwart, S., van Elteren, A., Whitehead, A., 2012\.  Simulations of Dense Stellar Systems with the AMUSE Software Toolkit, in: Capuzzo-Dolcetta, R., Limongi, M., Tornambè, A. (Eds.), Advances in Computational Astrophysics: Methods, Tools, and Outcome, p. 129.  arXiv:1111.3987.
  * McMillan (1986) McMillan, S.L.W., 1986.  The Vectorization of Small-N Integrators, in: Hut, P., McMillan, S.L.W. (Eds.), The Use of Supercomputers in Stellar Dynamics. volume 267, p. 156.  doi:.
  * Miller (1964) Miller, R.H., 1964.  Irreversibility in Small Stellar Dynamical Systems.  ApJ 140, 250.  doi:.
  * Miyamoto and Nagai (1975) Miyamoto, M., Nagai, R., 1975\.  Three-dimensional models for the distribution of mass in galaxies.  Publ. Astr. Soc. Japan 27, 533–543.
  * Navarro et al. (1996) Navarro, J.F., Frenk, C.S., White, S.D.M., 1996.  The Structure of Cold Dark Matter Halos.  ApJ 462, 563.  doi:, arXiv:astro-ph/9508025.
  * Pelupessy et al. (2012) Pelupessy, F.I., Jänes, J., Portegies Zwart, S., 2012.  N-body integrators with individual time steps from Hierarchical splitting.  New Astron. 17, 711–719.  doi:, arXiv:1205.5668.
  * Pelupessy and Portegies Zwart (2012) Pelupessy, F.I., Portegies Zwart, S., 2012\.  The evolution of embedded star clusters.  MNRAS 420, 1503–1517.  doi:, arXiv:1111.0992.
  * Pelupessy et al. (2013) Pelupessy, F.I., van Elteren, A., de Vries, N., McMillan, S.L.W., Drost, N., Portegies Zwart, S.F., 2013\.  The Astrophysical Multipurpose Software Environment.  A&A 557, A84.  doi:, arXiv:1307.3016.
  * Peters et al. (2018) Peters, B., Baniasadi, M., Baniasadi, M., Besseron, X., Estupinan Donoso, A., Mohseni, M., Pozzetti, G., 2018.  The XDEM Multi-physics and Multi-scale Simulation Technology: Review on DEM-CFD Coupling, Methodology and Engineering Applications.  arXiv e-prints , arXiv:1808.08028doi:, arXiv:1808.08028.
  * Portegies Zwart (2018) Portegies Zwart, S., 2018.  Computational astrophysics for the future.  Science 361, 979–980.  doi:, arXiv:1809.02600.
  * Portegies Zwart et al. (2009) Portegies Zwart, S., McMillan, S., Harfst, S., Groen, D., Fujii, M., Nualláin, B.Ó., Glebbeek, E., Heggie, D., Lombardi, J., Hut, P., Angelou, V., Banerjee, S., Belkus, H., Fragos, T., Fregeau, J., Gaburov, E., Izzard, R., Jurić, M., Justham, S., Sottoriva, A., Teuben, P., van Bever, J., Yaron, O., Zemp, M., 2009\.  A multiphysics and multiscale software environment for modeling astrophysical systems.  New Astron. 14, 369–378.  doi:, arXiv:0807.1996.
  * Portegies Zwart et al. (2013) Portegies Zwart, S., McMillan, S.L.W., van Elteren, E., Pelupessy, I., de Vries, N., 2013.  Multi-physics simulations using a hierarchical interchangeable software interface.  Computer Physics Communications 184, 456–468.  doi:, arXiv:1204.5522.
  * Portegies Zwart et al. (2020a) Portegies Zwart, S., Pelupessy, I., Martínez-Barbosa, C., van Elteren, A., McMillan, S., 2020a.  Non-intrusive hierarchical coupling strategies for multi-scale simulations in gravitational dynamics.  Communications in Nonlinear Science and Numerical Simulations 85, 105240.  doi:, arXiv:2002.11206.
  * Portegies Zwart et al. (2020b) Portegies Zwart, S., Pelupessy, I., Martínez-Barbosa, C., van Elteren, A., McMillan, S., 2020b.  Non-intrusive hierarchical coupling strategies for multi-scale simulations in gravitational dynamics.  Communications in Nonlinear Science and Numerical Simulations 85, 105240.  doi:, arXiv:2002.11206.
  * Portegies Zwart and Verbunt (1996) Portegies Zwart, S.F., Verbunt, F., 1996\.  Population synthesis of high-mass binaries.  A&A 309, 179–196.
  * Rin et al. (2017) Rin, R., Tomin, P., Garipov, T., Voskov, D., Tchelepi, H. (Eds.), 2017.  General Implicit Coupling Framework for Multi-Physics Problems. volume Day 3 Wed, February 22, 2017 of SPE Reservoir Simulation Conference.  URL: https://doi.org/10.2118/182714-MS, doi:, arXiv:https://onepetro.org/spersc/proceedings-pdf/17RSC/3-17RSC/D031S012R001/1298640/spe-182714-ms.pdf.
  * Strang (1968) Strang, G., 1968.  On the Construction and Comparison of Difference Schemes.  SIAM Journal on Numerical Analysis 5, 506–517.  doi:.
  * Terlevich (1987) Terlevich, E., 1987.  Evolution of N-body open clusters.  MNRAS 224, 193–225.  doi:.
  * Toonen et al. (2012) Toonen, S., Nelemans, G., Portegies Zwart, S., 2012.  Supernova Type Ia progenitors from merging double white dwarfs. Using a new population synthesis model.  A&A 546, A70.  doi:, arXiv:1208.6446.
  * van Elteren et al. (2019) van Elteren, A., Portegies Zwart, S., Pelupessy, I., Cai, M.X., McMillan, S.L.W., 2019.  Survivability of planetary systems in young and dense star clusters.  A&A 624, A120.  doi:, arXiv:1902.04652.
  * Veen and Hoekstra (2020) Veen, L.E., Hoekstra, A.G., 2020\.  Easing multiscale model design and coupling with muscle 3, in: Krzhizhanovskaya, V.V., Závodszky, G., Lees, M.H., Dongarra, J.J., Sloot, P.M.A., Brissos, S., Teixeira, J. (Eds.), Computational Science – ICCS 2020, Springer International Publishing, Cham. pp. 425–438.
  * Verlet (1967) Verlet, L., 1967.  Computer ”experiments” on classical fluids. i. thermodynamical properties of lennard-jones molecules.  Phys. Rev. 159, 98–103.  URL: http://link.aps.org/doi/10.1103/PhysRev.159.98, doi:.
  * Wielen (1967) Wielen, R., 1967.  Dynamical Evolution of Star Cluster Models, I.  Veroeffentlichungen des Astronomischen Rechen-Instituts Heidelberg 19, 1.
  * Xianmeng et al. (2021) Xianmeng, W., Mingyu, W., Xiao, H., Zhaoshun, W., Yinyu, C., Xu, L., Suxuan, G., 2021.  Multi-physics coupling simulation in virtual reactors.  Simulation 97, 687–702.  URL: https://doi.org/10.1177/0037549719881204, doi:.