# Working principles

This page provides an introduction to FONSim’s internal workings. Its aim is to provide the reader with an overview such that they can decide whether FONSim is the correct tool for their use case as well as reason about how they can use it the most effectively. Readers are also welcome to read the FONSim paper.

This document deliberately does not go into specifics or practicalities. For practical how-to guides on how to use FONSim, please have a look at the tutorials. There is also the formal code documentation which discusses all important functions and classes in FONSim. This page assumes that you have read the first paragraph of the introduction.

## Variables and terminals

A *system* in FONSim is defined as a bidirectional graph of *components*
connected by their *terminals*
and each terminal has one or more *variables*.
Two terminals are considered to be connected
when they are connected to a common *node*.

A *variable* corresponds with a scalar time dependent value,
for example a mass flow or a pressure.
Following classical systems modeling,
a variable in a terminal
is considered to be *over* or *through*.
Variables are also used to denote a component state
(for example, the amount of fluid inside).

A *terminal* is a set of variables
that together sufficiently describe
the interaction trough that terminal.
For example, for isothermal fluidics
a possible choice is a pressure (over variable)
and a mass flow (through variable).
As another example, mechanical rotational interactions
(such as a drive shaft)
can be described by torque (over variable)
and rotational velocity (through variable).

## Components and their equations

A *component* has one or more terminals as well as two functions.
The terminals define the component’s possibilities
to interact with other components,
while the functions describe its intrinsic behavior.
These two functions are the *update_state* function
(unless the component is stateless)
and the *evaluate* function.
They are expressed in plain Python
and can therefore contain any desired Python expression.

Both functions take as arguments the variable values.
The *update_state* function takes as additional arguments
the current state and the timestep size and it returns a new state.
(Any continuous differential equations
must therefore be discretized to a difference equation first.)
The *evaluate* function returns the equation residuals
that should equal zero at each simulation step.

Derivatives can be defined explicitly, or can be left out in which case FONSim will estimate them using finite differences. A more practical introduction to components that goes deeper into the usage specifics is provided by this tutorial.

## Simulation and solving

Once a *system* has been constructed,
FONSim’s *simulation* module allows
to derive numerical results from it.
The default FONSim solver is a fairly straightforward implementation
of Newton-Raphson
with backward Euler discretization for improved stability.
(Given that this backward discretization
is not apparent in the component definition,
it is relatively easy to employ forward Euler instead.)
In addition, the default solver is equipped with automatic timestep adjustment.
In usage, simulating a defined system takes minimal effort
(see, for example, this tutorial)
because FONSim makes abstraction of the simulation procedure.

### Initialization

When asked to simulate a system, FONSim initializes the candidate solution vector \(\mathbf{x_a}\), which contains the candidate values of all variables in terminals, and the square system matrix \(\mathbf{J}\). Hashtables (Python dictionaries) map variables in the object oriented system structure to matrix indices in \(O(1)\) time. FONSim also builds using Kirchoff laws the network equations that describe the relations between variables of different components. These equations are linear and constant for a given system.

The total residual vector \(\mathbf{f}\) takes as argument the candidate solution vector \(\mathbf{x}\) and returns the residual. If the solution were to be exactly correct, the residual would be a vector of zeros. The system matrix \(\mathbf{J}\) is the derivative of this residuals vector \(\mathbf{f}\) to the solution vector \(\mathbf{x}\) (i.e. Jacobian matrix).

The upper part \(\mathbf{J_n}\) of the system matrix expresses the network equations and is therefore constant over a simulation. On the other hand, the lower part \(\mathbf{J_a}\) expresses the derivatives of the component equations. As those equations tend to be nonlinear, \(\mathbf{J_a}\) has to be updated in every iteration of the Newton outer loop.

### Run

After these preparations,
FONSim starts iterating over the timesteps,
solving timestep by timestep.
At each timestep FONSim solves
the nonlinear residual function \(\mathbf{f}\)
using the Newton-Raphson method.
Each Newton outer loop iteration
starts with evaluating this residual vector
with the candidate solution \(\mathbf{x}\)
For the component residuals, this is done by iterating over all components
and evaluating their *update_state* and *evaluate* functions.
(FONSim does not do any symbolic math on the given expressions.)
The network residuals follow from multiplying
\(\mathbf{J_n}\) with \(\mathbf{x}\).
Subsequently, this residual is backprojected onto \(\mathbf{x}\)
using a highly optimized Numpy linear system solver.
This outer loop is repeated until all residuals
are lower than set thresholds.
(These thresholds depend on the variables’ ranges
to avoid scale effects.)

While evaluating the component residuals of each component,
FONSim also gathers the derivatives
of the *evaluate* and the *update_state* functions.
The product of these two is then written
to the appropriate location in \(\mathbf{J_a}\),
Had forward Euler been used, only the derivatives
of the *evaluate* function would be used.

The timestep adjustment is handled using an interpolation criterion. Given \(\mathbf{x}\) solved at three subsequent simulation steps \(k-2\), \(k-1\) and \(k\), the vector \(\mathbf{x}^{k-1}\) must be within the desired precision from the \((\mathbf{x}^{k-2} + \mathbf{x}^{k})/2\). If this is not the case, the timestep is reduced, while if it is far closer than the desired precision, the timestep is increased. This desired precision can be specified as relative to the variable’s range. Each time the timestep is reduced, the simulator returns two steps to ensure proper sampling.

Once the simulation has finished, the simulation results are easily accessed for visualization or other purposes. One can access them via the components because their internal variables are mapped to the simulation solution matrix.

## Conclusion

This page explains the working principles of FONSim. The behavior of the interconnected components is defined using variables, terminals and equations. Together the connected components form a system, which is modeled as a bidirectional graph. Components have terminals that define how they can interact with other components and functions that describe their intrinsic behavior, which can be expressed in plain Python. FONSim’s default solver is a straightforward implementation of the Newton-Raphson method with backward Euler discretization, and it has automatic timestep adjustment. The system matrix expresses the network equations and the derivatives of component equations. Simulating a defined system takes minimal effort because FONSim makes abstraction of the simulation procedure.

For further reading, it is suggested to have a look at the tutorials and the formal code documentation.