Finesse internal code structure¶

Interferometer Matrix Equation¶

(21)¶\[\biggl( \begin{matrix} \text{interferometer} \\ \text{matrix} \end{matrix} \biggl) \times \biggl( \vec{x}_{\rm sol} \biggl) = \biggl( \vec{x}_{\rm RHS} \biggl)\]

The whole of Finesse (excluding things like the parser, and beamtracing/propagation code) revolves around setting up and solving the matrix equation (21) by inverting the interferometer matrix and multiplying it with \(\vec{x}_{\rm RHS}\)

There is a separate matrix for the carrier and signal calculations

Most of this page mentions code currently implementing the sparse solver, which use Suitesparse KLU sparse matrix solver under the hood. In principle the code architecture supports implementing different solvers, like a possible graph-based symbolic solver.

Solution Vector¶

The desired output is the solution vector \(\vec{x}_{\rm sol}\). In the codebase this vector is named out_view in simulations/homsolver.pxd. It contains the field (light amplitude and phase) at every Node, for every frequency and for every higher order mode.

The index into this vector is calculated in the finesse.simulations.homsolver.HOMSolver.field() method in simulations/homsolver.pyx, which is used in the finesse.simulations.homsolver.HOMSolver.get_out() to get a single value and finesse.simulations.homsolver.HOMSolver.node_field_vector() to get a pointer to the HOM vector.

Most of the detectors end up calling one of these functions (indirectly, via their workspaces) to get the field at a certain node and do some further calculations. See for instance c_ad_single_field_output in detectors/compute/amplitude.pyx.

Right hand side vector¶

The right hand side vector \(\vec{x}_{\rm RHS}\) has the same structure as the solution vector, but all the values are initially set to zero, except for source nodes such as lasers and modulators (described in more detail in Input fields or the ‘right hand side’ vector).

Source values are set using the finesse.simulations.sparse.solver.SparseSolver.set_source() and set_source_fast functions defined in simulations/sparse/solver.pyx. You can see it is being indexed with the same field and field_fast functions mentioned above.

You can see an example in c_laser_carrier_fill_rhs in components/modal/laser.pyx, where the laser field values are being set.

Interferometer Matrix¶

The interferometer matrix contains the linear equations that couple the field values of different nodes together.

The following simple example describes the field couplings for a mirror (also described here Mirror):

(22)¶\[\left( \begin{array}{c} \rm Out1 \\ \rm Out2 \end{array} \right)=\left( \begin{array}{cc} a_{11} &a_{21} \\ a_{12} & a_{22} \end{array}\right) \left( \begin{array}{c} \rm In1 \\ \rm In2 \end{array}\right)\]

(23)¶\[\left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ -a_{11} & 1 & -a_{21} & 0 \\ 0 & 0 & 1 & 0 \\ -a_{12} & 0 &-a_{22} & 1 \end{array}\right)\times \left(\begin{array}{c} {\rm In1} \\ {\rm Out1} \\ {\rm In2}\\ {\rm Out2}\\\end{array}\right)= \left(\begin{array}{c} {\rm In1} \\ 0 \\ {\rm In2}\\ 0\\\end{array}\right).\]

In a typical interferometer, most components only couple with components directly attached to them. Therefore the interferometer matrix will be ‘sparse’ (most entries are zero).

CCSMatrix¶

For memory efficiency and performance, finesse stores the IFO matrix in compressed column storage format (CCS/CSC), which is implemented in the finesse.cymath.cmatrix.CCSMatrix class in cymath/cmatrix.pyx. This resource has a nice explanation of the format.

The most important methods are:

finesse.cymath.cmatrix.CCSMatrix.declare_equations()

finesse.cymath.cmatrix.CCSMatrix.declare_submatrix_view()

finesse.cymath.cmatrix.CCSMatrix.construct()

Multiple calls to declare_equations define the structure and the order of the sparse matrix. Since the sparse matrix represents a linear system of equations, every row relates to an ‘equation’ coupling the corresponding field to other fields in the system of equations, with the values in the row describing the coupling coefficients. Equation (23) describes the following 4 equations:

\[\begin{aligned} & \mathrm{In1} && = \mathrm{In1} \\ -a_{11} \cdot & \mathrm{In1} + \mathrm{Out1} - a_{21} \cdot \mathrm{In2} && = 0 \\ & \mathrm{In2} && = \mathrm{In2} \\ -a_{12} \cdot & \mathrm{In1} + \mathrm{Out2} - a_{22} \cdot \mathrm{In2} && = 0 \\ \end{aligned}\]

We can set up this problem using finesse.cymath.cmatrix.KLUMatrix (which derives from finesse.cymath.cmatrix.CCSMatrix, but interfaces to KLU as well so we can solve the matrix at the end):

import numpy as np
from finesse.cymath.cmatrix import KLUMatrix

mat = KLUMatrix(name="mirror_example")
In1 = mat.declare_equations(Neqs=1, index=0, name="I_In1", is_diagonal=True)
Out1 = mat.declare_equations(Neqs=1, index=1, name="I_Out1", is_diagonal=True)
In2 = mat.declare_equations(Neqs=1, index=2, name="I_In2", is_diagonal=True)
Out2 = mat.declare_equations(Neqs=1, index=3, name="I_Out2", is_diagonal=True)

mat.construct()
mat.print_matrix()

Matrix mirror_example: nnz=4 neqs=4
    (col, row) = value
    (0, 0) = (1+0j) : I_In1 mode=0 -> I_In1 mode=0
    (1, 1) = (1+0j) : I_Out1 mode=0 -> I_Out1 mode=0
    (2, 2) = (1+0j) : I_In2 mode=0 -> I_In2 mode=0
    (3, 3) = (1+0j) : I_Out2 mode=0 -> I_Out2 mode=0

Here we also see the function of finesse.cymath.cmatrix.CCSMatrix.construct(), which needs to be called after declaring the equations to allocate the memory and fill in the values.

So far we have only set up a system of 4 equations where every field only directly couples to itself (setting is_diagonal=True fills in a 1 in the diagonal positions). Note that by setting Neqs=1 we only declare a single equations per index, but we can also declare a block of equations (representing higher order modes for example).

To set up coupling between the different fields, we use finesse.cymath.cmatrix.CCSMatrix.declare_submatrix_view():

Note

To make the matrix printing work correctly, we have to comply with the convention that the names of submatrices created by finesse.cymath.cmatrix.CCSMatrix.declare_equations() start with an I, while submatrices created by finesse.cymath.cmatrix.CCSMatrix.declare_submatrix_view() should not start with an I.

a11 = mat.declare_submatrix_view(In1.from_idx, Out1.from_idx, "_In1->Out1", False)
a22 = mat.declare_submatrix_view(In2.from_idx, Out2.from_idx, "_In2->Out2", False)

a12 = mat.declare_submatrix_view(In1.from_idx, Out2.from_idx, "_In1->Out2", False)
a21 = mat.declare_submatrix_view(In2.from_idx, Out1.from_idx, "_In2->Out1", False)

mat.construct()
mat.print_matrix()

print(f"{a11.shape=}")

Matrix mirror_example: nnz=8 neqs=4
    (col, row) = value
    (0, 0) = (1+0j) : I_In1 mode=0 -> I_In1 mode=0
    (0, 1) = 0j     : I_In1 mode=0 -> I_Out1 mode=0
    (0, 3) = 0j     : I_In1 mode=0 -> I_Out2 mode=0
    (1, 1) = (1+0j) : I_Out1 mode=0 -> I_Out1 mode=0
    (2, 1) = 0j     : I_In2 mode=0 -> I_Out1 mode=0
    (2, 2) = (1+0j) : I_In2 mode=0 -> I_In2 mode=0
    (2, 3) = 0j     : I_In2 mode=0 -> I_Out2 mode=0
    (3, 3) = (1+0j) : I_Out2 mode=0 -> I_Out2 mode=0

a11.shape=(1, 1)

The first two arguments are the indices of the nodes between which we want to declare a coupling. In this case, since there is only 1 equations per node, we are creating a memory view to a 1x1 matrix, which represents a single coupling coefficient. The power of this abstraction is more apparent when you have multiple equations per node, for which declare_submatrix_view returns a view of a dense matrix of which the elements are located at the appropriate locations in the overall sparse matrix.

We can now ‘fill’ the sparse matrix by overwriting the values in the submatrix views using the NumPy array slicing syntax:

r = np.sqrt(0.9)
t = np.sqrt(0.1)

a11[:] = r
a22[:] = r
a12[:] = 1j * t
a21[:] = 1j * t

mat.print_matrix()

# print in dense format
np.set_printoptions(precision=1, linewidth=250)
print(mat.to_scipy_csc().todense())

Matrix mirror_example: nnz=8 neqs=4
    (col, row) = value
    (0, 0) = (1+0j)                    : I_In1 mode=0 -> I_In1 mode=0
    (0, 1) = (-0.9486832980505138+0j)  : I_In1 mode=0 -> I_Out1 mode=0
    (0, 3) = (-0-0.31622776601683794j) : I_In1 mode=0 -> I_Out2 mode=0
    (1, 1) = (1+0j)                    : I_Out1 mode=0 -> I_Out1 mode=0
    (2, 1) = (-0-0.31622776601683794j) : I_In2 mode=0 -> I_Out1 mode=0
    (2, 2) = (1+0j)                    : I_In2 mode=0 -> I_In2 mode=0
    (2, 3) = (-0.9486832980505138+0j)  : I_In2 mode=0 -> I_Out2 mode=0
    (3, 3) = (1+0j)                    : I_Out2 mode=0 -> I_Out2 mode=0

[[ 1. +0.j   0. +0.j   0. +0.j   0. +0.j ]
 [-0.9+0.j   1. +0.j   0. -0.3j  0. +0.j ]
 [ 0. +0.j   0. +0.j   1. +0.j   0. +0.j ]
 [ 0. -0.3j  0. +0.j  -0.9+0.j   1. +0.j ]]

Note

The submatrix views deriving from finesse.cymath.cmatrix.SubCCSView also provide optimised c-functions for filling values. The z in the naming refers to complex values, a to a scalar and M to a matrix.

fill_za fills the entire submatrix with the same value a
fill_zm copies over values from matrix M with the same shape
fill_za_zm fills the submatrix with the scalar-matrix product of a * M

The reason they all exist as separate functions is for optimization only.

These submatrix view objects are being kept track of in ‘connections’ objects like finesse.components.modal.mirror.MirrorOpticalConnections, which are being stored in Workspace (like finesse.components.modal.mirror.MirrorWorkspace)so eventually the workspaces can call their fill functions to fill in couplings in the matrix. See for example the mirror_fill_optical_2_optical in components/modal/mirror.pyx, which ends up calling the fill_negative_za_zm_2 method of SubCCSView.

We have now created a matrix identical to (23), and we can get to ‘solving’ it. This consists of three steps:

Factoring the matrix into LU form with finesse.cymath.cmatrix.KLUMatrix.factor()
Inverting the matrix with finesse.cymath.cmatrix.KLUMatrix.solve()
Multiplying the inverse matrix with \(\vec{x}_{\rm RHS}\) to get \(\vec{x}_{\rm sol}\). This is also done by finesse.cymath.cmatrix.KLUMatrix.solve(), which returns a view of \(\vec{x}_{\rm RHS}\) which has been overwritten with the solution values

We will first inject an initial value for the incoming field:

mat.set_rhs(In1.from_idx, 1.0)

And finally do the solving:

mat.print_rhs()

mat.factor()
mat.solve()

mat.print_rhs()

Vector mirror_example: neqs=4
    (row) = value
    (0) = (1+0j) : I_In1 mode=0
    (1) = 0j     : I_Out1 mode=0
    (2) = 0j     : I_In2 mode=0
    (3) = 0j     : I_Out2 mode=0


Vector mirror_example: neqs=4
    (row) = value
    (0) = (1+0j)                  : I_In1 mode=0
    (1) = (0.9486832980505138+0j) : I_Out1 mode=0
    (2) = 0j                      : I_In2 mode=0
    (3) = 0.31622776601683794j    : I_Out2 mode=0

You can see the matrix in action by passing debug_mode=True via the simulation_options dict of finesse.model.Model.run

import finesse
import numpy as np

model = finesse.Model()

np.set_printoptions(linewidth=200, precision=1)

model.parse(
    """
        l l1
        s s1 l1.p1 m1.p1
        m m1 R=0.9 T=0.1
    """
)
model.run(simulation_options={"debug_mode": True});

Matrix carrier: nnz=12 neqs=6
    (col, row) = value
    (0, 0) = (1+0j)                   : I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0 -> I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (1, 1) = (1+0j)                   : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0 -> I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0
    (1, 2) = (-1-0j)                  : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0 -> I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (2, 2) = (1+0j)                   : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (2, 3) = (-0.9486832980505138-0j) : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (2, 5) = -0.31622776601683794j    : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0
    (3, 0) = (-1-0j)                  : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0 -> I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (3, 3) = (1+0j)                   : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4, 3) = -0.31622776601683794j    : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4, 4) = (1+0j)                   : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0
    (4, 5) = (-0.9486832980505138+0j) : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0
    (5, 5) = (1+0j)                   : I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0


Matrix in dense format:

[[ 1. +0.j   0. +0.j   0. +0.j  -1. +0.j   0. +0.j   0. +0.j ]
 [ 0. +0.j   1. +0.j   0. +0.j   0. +0.j   0. +0.j   0. +0.j ]
 [ 0. +0.j  -1. +0.j   1. +0.j   0. +0.j   0. +0.j   0. +0.j ]
 [ 0. +0.j   0. +0.j  -0.9+0.j   1. +0.j   0. -0.3j  0. +0.j ]
 [ 0. +0.j   0. +0.j   0. +0.j   0. +0.j   1. +0.j   0. +0.j ]
 [ 0. +0.j   0. +0.j   0. -0.3j  0. +0.j  -0.9+0.j   1. +0.j ]]

Right hand side vector:

Vector carrier: neqs=6
    (row) = value
    (0) = (1.3416407864998738+0j) : I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (1) = (1.4142135623730951+0j) : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0
    (2) = (1.4142135623730951+0j) : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (3) = (1.3416407864998738+0j) : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4) = 0j                      : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0
    (5) = 0.447213595499958j      : I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0

Note that you can have direct access to the finesse.simulations.sparse.simulation.SparseMatrixSimulation instance by using the context manager finesse.model.Model.built(), and thereby to the finesse.simulations.sparse.KLU.KLUSolver and eventually the finesse.cymath.cmatrix.KLUMatrix instance:

with model.built() as sim:
    carrier_solver = sim.carrier
    klu_mat = carrier_solver.M()
    klu_mat.print_matrix()
    klu_mat.print_rhs()

Matrix carrier: nnz=12 neqs=6
    (col, row) = value
    (0, 0) = (1+0j)                   : I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0 -> I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (1, 1) = (1+0j)                   : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0 -> I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0
    (1, 2) = (-1-0j)                  : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0 -> I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (2, 2) = (1+0j)                   : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (2, 3) = (-0.9486832980505138-0j) : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (2, 5) = -0.31622776601683794j    : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0
    (3, 0) = (-1-0j)                  : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0 -> I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (3, 3) = (1+0j)                   : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4, 3) = -0.31622776601683794j    : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4, 4) = (1+0j)                   : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0
    (4, 5) = (-0.9486832980505138+0j) : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0
    (5, 5) = (1+0j)                   : I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0 -> I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0


Vector carrier: neqs=6
    (row) = value
    (0) = (1.3416407864998738+0j) : I,node=l1.p1.i,f=0,fidx=0,Neq=1 mode=0
    (1) = (1.4142135623730951+0j) : I,node=l1.p1.o,f=0,fidx=1,Neq=1 mode=0
    (2) = (1.4142135623730951+0j) : I,node=m1.p1.i,f=0,fidx=2,Neq=1 mode=0
    (3) = (1.3416407864998738+0j) : I,node=m1.p1.o,f=0,fidx=3,Neq=1 mode=0
    (4) = 0j                      : I,node=m1.p2.i,f=0,fidx=4,Neq=1 mode=0
    (5) = 0.447213595499958j      : I,node=m1.p2.o,f=0,fidx=5,Neq=1 mode=0

KLUMatrix¶

As seen above, the KLUMatrix class, implements methods to factor and solve sparse matrices by linking into KLU. KLU is a library included in SuiteSparse which stand for “Clark Kent” LU factorization algorithm (as opposed to “SuperLU”).

This is also the source of the infamous “KLU_singular error”, which can easily be replicated by passing a singular matrix:

from finesse.cymath.cmatrix import KLUMatrix

mat = KLUMatrix(name="singular_example")
In1 = mat.declare_equations(Neqs=1, index=0, name="I_In1", is_diagonal=True)
Out1 = mat.declare_equations(Neqs=1, index=1, name="I_Out1", is_diagonal=True)
a12 = mat.declare_submatrix_view(In1.from_idx, Out1.from_idx, "_In1->Out1", False)
a21 = mat.declare_submatrix_view(Out1.from_idx, In1.from_idx, "_Out1->In1", False)

mat.construct()

a12[:] = 1
a21[:] = 1

mat.factor()

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[8], line 14
     11 a12[:] = 1
     12 a21[:] = 1
---> 14 mat.factor()

File src/finesse/cymath/cmatrix.pyx:1736, in finesse.cymath.cmatrix.KLUMatrix.factor()

File src/finesse/cymath/cmatrix.pyx:1771, in finesse.cymath.cmatrix.KLUMatrix.factor()

Exception: An error occurred in KLU during factoring: STATUS=1 (KLU_SINGULAR)

Matrix singular_example: nnz=4 neqs=2
    (col, row) = value
    (0, 0) = (1+0j)  : I_In1 mode=0 -> I_In1 mode=0
    (0, 1) = (-1+0j) : I_In1 mode=0 -> I_Out1 mode=0
    (1, 0) = (-1+0j) : I_Out1 mode=0 -> I_In1 mode=0
    (1, 1) = (1+0j)  : I_Out1 mode=0 -> I_Out1 mode=0

As you can see, at the core the matrix solving is rather straightforward. Most of the complexity of Finesse arises from keeping track of all the state: the components in the model, their changing parameters, the different frequencies and higher order modes and the connections between the components.

Simulation¶

The simulation class takes the Model object and ‘builds’ a simulation from it. This includes the following:

Beam Tracing functionality (this can be done without solving the matrix equation)
Storing common simulation settings
Checks is a signal simulation is required
Checks which Parameters are changing and update symbolic expressions
Determines is the simulation is modal
Keep track of the detector workspaces

The BaseSimulation class is found in simulations/simulation.pyx and the child class SparseMatrixSimulation in simulations/sparse/simulation.pyx. The separation is such that Finesse might support different types of solvers (like a graph-based solver) in the feature.

Most importantly, it builds the carrier Solver and (if necessary) the signal Solver objects. The simulation object is saved in the finesse.analysis.actions.base.AnalysisState object, which in turn gets passed to the finesse.analysis.actions.base.Action._do method of the actions.

The actions then usually call state.sim.run_carrier or state.sim.run_signal, depending on what the action needs for it’s output.

You can manually create a SparseMatrixSimulation object using the finesse.model.Model.built() context manager:

Warning

Manually running simulations like this is only for demonstration the internals of Finesse!

import numpy as np
import finesse


model = finesse.Model()

model.parse(
    """
    l l1 P=1
    s s1 l1.p1 m1.p1
    m m1 R=0.1 T=0.9
    """
)
with model.built() as sim:
    sim.run_carrier()
    for node, val in zip(
        sim.carrier.nodes.keys(), np.array(sim.carrier.M().rhs_view[0, :]),
        strict=True
    ):
        print(node, val)

l1.p1.i (0.447213595499958+0j)
l1.p1.o (1.4142135623730951+0j)
m1.p1.i (1.4142135623730951+0j)
m1.p1.o (0.447213595499958+0j)
m1.p2.i 0j
m1.p2.o 1.3416407864998738j

This example is loosely based on the finesse.analysis.actions.dc.DCFields class.

Solver¶

The solver object couples the nodes from the Model to the underlying matrix solver. In principle it would support different kinds of matrix solvers, but currently only the KLUMatrix class is used.

The inheritance tree is BaseSolver->HOMSolver->SparseSolver->KLUSolver. The classes are implemented in simulations and simulations/sparse.

ModelElement¶

Detector¶

Port¶

Node¶

Parameters¶

Workspace¶

Model¶

Signal Matrix¶

Click to download example as python script

Click to download example as Jupyter notebook