Class: Base system
This class simply initiates a system, defines the grid and contains the basic functionality for evolving in time.
Types of functions
There are five different types of functions:
conf_
functions: Configures the state of the system, for instance by setting an initial condition. Output nothing.evolve_
functions: Evolves the state in time according to some equation of motion.calc_
functions: Calculates something from the state, returns that which has been calculated.plot_
functions: Functions tailored to plot specific things. Output the axes and figure.get_
functions: Functions that return a component of a tensor field. Relevant for symmetric and antisymmetric tensors where it is not convenient to save all elements.
General keywords and parameters
The only required input argument to the BaseSystem class is the dim
argument, which specifies the dimension of the system.
In some cases, default values of other parameter depend on the value of dim
, and are represented by curly brackets:
These are the optional keywords for the BaseSystem
class.
Keyword  Definition  Default value 

xmin 
Minimum value of \(x\) of the simulation domain  \(0\) 
ymin 
Minimum value of \(y\) of the simulation domain  \(0\) 
zmin 
Minimum value of \(z\) of the simulation domain  \(0\) 
xmax 
Maximum value of \(x\) of the simulation domain.  \(100\) 
ymax 
Maximum value of \(y\) of the simulation domain  \(\left \lbrace \begin{array}{c} 1 \\ 100 \\ 100 \\ \end{array} \right \rbrace\) 
zmax 
Maximum value of \(z\) of the simulation domain  \(\left \lbrace \begin{array}{c} 1 \\ 1 \\ 100 \\ \end{array} \right \rbrace\) 
xRes 
Resolution of the \(x\) axis  \(101\) 
yRes 
Resolution of the \(y\) axis  \(\left \lbrace \begin{array}{c} 1 \\ 101 \\ 101 \\ \end{array} \right \rbrace\) 
zRes 
Resolution of the \(z\) axis  \(\left \lbrace \begin{array}{c} 1 \\ 1 \\ 101 \\ \end{array} \right \rbrace\) 
dx 
Spacing between points on the \(x\)axis. Trumps xRes if provided. xmax will be modified to match. 
\(\frac{\texttt{xmax}\texttt{xmin}}{\texttt{xRes}} = 1\) 
dy 
Spacing between points on the \(y\)axis. Trumps yRes if provided. 
\(\frac{\texttt{ymax}\texttt{ymin}}{\texttt{yRes}} = 1\) 
dz 
Spacing between points on the \(x\)axis. Trumps zRes if provided. 
\(\frac{\texttt{zmax}\texttt{zmin}}{\texttt{zRes}} = 1\) 
xlim 
List or tuple consisting of the lower and upper limit for the simulation domain in the \(x\)direction. Trumps xmin and xmax if provided. 
\((\texttt{xmin},\texttt{xmax}) = (0,101)\) 
ylim 
List or tuple consisting of the lower and upper limit for the simulation domain in the \(y\)direction. Trumps ymin and ymax if provided. 
\((\texttt{ymin},\texttt{ymax}) = \left \lbrace \begin{array}{c} (0,1) \\ (0,101) \\ (0,101) \\ \end{array} \right \rbrace\) 
zlim 
List or tuple consisting of the lower and upper limit for the simulation domain in the \(z\)direction. Trumps zmin and zmax if provided. 
\((\texttt{xmin},\texttt{xmax}) = \left \lbrace \begin{array}{c} (0,1) \\ (0,1) \\ (0,101) \\ \end{array} \right \rbrace\) 
time 
Float specifying the time of initialization  \(0\) 
a0 
Characteristic length scale associated with the system, in units of which all plots will be scaled. This is also the default width used with the coarsegraining operation.  \(1\) 
X 
2D numpy array with position coordinates. Typically used when the coordinates are not a regular grid.  None 
Y 
2D numpy array with position coordinates. Typically used when the coordinates are not a regular grid.  None 
Z 
2D numpy array with position coordinates. Typically used when the coordinates are not a regular grid.  None 
From these keywords, a number of useful parameters are constructed, given in the table below.
Parameter  Definition  Value 

x 
Numpy array with dimensions \(\left \lbrace \begin{array}{l} \texttt{xRes} \\ \texttt{xRes}\times 1 \\ \texttt{xRes}\times 1 \times 1 \\ \end{array} \right \rbrace\) consisting of the grid points from xmin to (including) xmaxdx . 

y 
Numpy array with dimensions \(\left \lbrace \begin{array}{l} 1 \\ 1 \times \texttt{yRes} \\ 1 \times \texttt{yRes} \times 1 \\ \end{array} \right \rbrace\) consisting of the grid points from ymin to (including) ymaxdy . 

z 
Numpy array with dimensions \(\left \lbrace \begin{array}{l} 1 \\ 1 \\ 1 \times 1 \times \texttt{zRes} \\ \end{array} \right \rbrace\) consisting of the grid points from zmin to (including) zmaxdz . 

xmidi 
Index of the mid \(x\)value. In the case of an odd xRes , this midpoint index will not hit the middle exactly but undershoot by dx/2 . 

xmid 
The \(x\) value given by xmidi . 

ymidi 
Index of the mid \(y\)value. In the case of an odd yRes , this midpoint index will not hit the middle exactly but undershoot by dy/2 . 

ymid 
The \(y\) value given by ymidi . 

zmidi 
Index of the mid \(z\)value. In the case of an odd zRes , this midpoint index will not hit the middle exactly but undershoot by dz/2 . 

zmid 
The \(z\) value given by zmidi . 

size_x 
Size of the \(x\)axis  \(\texttt{xmax}  \texttt{xmin}\) 
size_y 
Size of the \(y\)axis  \(\texttt{ymax}  \texttt{ymin}\) (1 if dim \(<2\)) 
size_z 
Size of the \(z\)axis  \(\texttt{zmax}  \texttt{zmin}\) (1 if dim \(<3\)) 
size_min 
Minimum value of the simulation domain  \(\left \lbrace \begin{array}{c} \texttt{size_x} \\ \texttt{min}(\texttt{size_x}, \texttt{size_y})\\ \texttt{min}(\texttt{size_x}, \texttt{size_y}, \texttt{size_z}) \\ \end{array} \right \rbrace\) 
size_max 
Maximum value of the simulation domain  \(\left \lbrace \begin{array}{c} \texttt{size_x} \\ \texttt{max}(\texttt{size_x}, \texttt{size_y})\\ \texttt{max}(\texttt{size_x}, \texttt{size_y}, \texttt{size_z}) \\ \end{array} \right \rbrace\) 
dV 
Volume element of the grid.  \(\left \lbrace \begin{array}{l} \texttt{dx} \\ \texttt{dx} \times \texttt{dy} \\ \texttt{dx} \times \texttt{dy} \times \texttt{dz} \\ \end{array} \right \rbrace\). 
volume 
Volume of the simulation domain  \(\texttt{size_x} \times \texttt{size_y} \times \texttt{size_z}\) 
Note that even though variables like yRes
, zRes
etc. are defined in cases where they are not relevant, such as for a \(1\)dimensional system, they play no significant role in any calculations in such situations.
Periodic boundary conditions means that xmax
and xmin
are identified as the same point.
Coarsegraining
A common and useful method is that of coarsegraining, which is defined as
where \(\mathcal K(\mathbf r'\mathbf r)\) is a Gaussian kernel given by
From a numerical point of view, this is done in Fourier space since, by the convolution theorem,
Thus, we need the Fourier transform of \(\mathcal K\), which is
This is why we have the following function
which calculates \(\mathcal K_{\mathfrak f}\).
Typically, a field is coarsegrained with a width using the following piece of code
The Gaussian function is actually so useful that is given by can be calculated using
Vortex fields
A general feature that will be reused is that of vortex fields. An angle field is a field where each point in space corresponds to an angle \(\theta \in \mathcal S^n\). A vortex is a topological defect in an angle field, around which the circulation is some integer multiple of the covering of \(\mathcal S^n\).
Angle field of a single vortex in two dimensions
In two dimensions, the angle field takes values \(\theta \in [\pi,\pi \rangle\) and a vortex is a point \(\mathbf r_0\). The angle field produced by the vortex has a circulation which is a multiple integer of \(2\pi\), i.e.,
where \(s_n\) is the integer charge of the vortex. A possible angle field for a vortex positioned at \((x_0,y_0)\) is given by
Angle field of a vortex ring in three dimensions
For a ring vortex in three dimensions centered at \(\mathbf{r_0}\) with radius \(R\), an angle field with the correct topological charge is given by first calculating the auxiliary quantities
and then calculating
These expressions are based on the geometry depicted in the following figure.
Vortex ring angle field explanation: Geometry of a vortex ring in the plane given by \(\vec n\). \(\mathcal N'\) is the plane normal to the tangent vector \(\vec t'\) at \(\vec r'\) upon which we impose a Cartesian coordinate system to determine the angles \(\theta_1\), \(\theta_2\) that are used to construct the (inset) initial angle field. Figure reprinted from Ref.^{3} with permission.
The angle field is then given by
and is implemented in the function calc_angle_field_vortex_ring
.
Periodic boundary conditions: Numerical implementation of angle fields
Apart from the angle field of a single vortex, the other fields are compatible with periodic boundary conditions. The expressions for these fields, however, are really only valid for an infinite region. When this is imposed on periodic boundary conditions, it results in spurious boundary effects, especially if either of the vortices is placed near the edge of the simulation domain. By simply inserting the vortices directly, we get what is shown in the following figure (a).
Numerical implementaiton of periodic angle fields: The angle field of panel (a) has been filtered by the field \(F\) with \(w=0.2x_{\textrm{max}}\) to produce the periodic field given in panel (c). This field is simply rolled to produce a different position for the dipole in panel (d).
This field is not periodic on the domain. This typically causes the unintentional nucleation of vortices and strain on the boundary. We therefore seek to modify the fields so that they don't "see" the periodic boundary conditions.
In order to produce a field that is periodic on the domain, we transform the field \(\theta\) to a complex field \(\eta = e^{i \theta}\). The argument of this complex field has the correct winding configuration of the vortex dipole. However, we want to smooth this field so that it goes to \(1\) (\(\theta=0)\) at the borders. To do so, we introduce the filter function
where \(r^2 = (xx_{\textrm{mid}})^2 + (yy_{\textrm{mid}})^2\), which is a function that is zero in the center region and goes to \(1\) at infinity over a width of \(w\). The filtered field \(\eta\) is then obtained by making a smooth function that goes from the previous angle field according to
\(\tilde \eta\) is shown in Figure (c). The value of \(w\) and \(R\) can be adjusted and are found in the source code as width
and radius
, respectively.
From this section, it is clear that the initial vortex dipole should not be too large. Thus, we have included a warning in case it is attempted to initiate a dipole with a larger distance than a certain threshold.
Numerical integration scheme
The systems of interest for this code are those that can be written on the form
where \(\omega\) is a linear differential operator and \(N\) is a nonlinear operator (function of \(\psi\)). The following table shows some examples from the models that we will discuss in the following chapters.
Model  \(\omega\)  \(\omega_{\mathfrak f}(\mathbf{k})\)  \(N\) 

Quantum Mechanics  $\frac{1}{2}i \nabla^2 $  \(\frac{1}{2} i \mathbf{k}^2\)  \( i V\) 
BEC  \((i+\gamma) (1+\frac{1}{2}\nabla^2)\)  \((i+\gamma) (1\frac{1}{2}\mathbf k^2)\)  \( (i + \gamma) (V_{ext} + \psi \psi^*)\psi\) 
Active Nematic  \(\frac{K}{\gamma} \nabla^2 +\frac{AB}{\gamma}\)  \(\frac{K}{\gamma} k^2 +\frac{AB}{\gamma}\)  \( \mathbf{u}\cdot \nabla Q + Q \Omega \Omega Q  \frac{2A}{\gamma}Q^2_{kk}Q\) 
Table: Examples of time evolution operators, nonconserved.
In the following, we will explain the method of evolution of exponential time differencing for stiff systems^{1}. This will result in two integration schemes, the exponential time differencing second order Runge Kutta 2 (ETD2RK) scheme and the forth order ETD4RK scheme. As in Ref. [coxExponentialTimeDifferencing2002], we will show an intuitive way to obtain the former and only recite the expressions for the latter.
The ETD2RK scheme
To see how we obtain the ETD2RK scheme, we take the Fourier transformation of Eq. ref:timeevolution
, and get:
Where \(\psi_{\mathfrak f}(\mathbf{k})\) is the Fourier transform of \(\psi(\mathbf{r})\). Integrating from \(t\) to \(t+ \Delta t\), one finds:
we multiply by \(e^{ \omega_{\mathfrak f}(t+ \Delta t)}\) and get:
This is an exact result, however, the last integral is unknown. In order to calculate the last integral here, we approximate it by \(N (t+\tau) \approx N_{\mathfrak f 0} + \frac{\Delta N_{\mathfrak f}}{\Delta t} \tau\) where \(N_{\mathfrak f 0} = (N(\psi(t))_{\mathfrak f}\) and \(\Delta N_{\mathfrak f} = N_{\mathfrak f}(t+\Delta t)N_{\mathfrak f}(t)\). We also change the integration limits from \(\tau \in [t,t+\Delta t]\) to \(\tau \in [0,\Delta t]\), which gives:
To find \(\psi_{\mathfrak f} (t+\Delta t)\), we would need to know the value \(N_{\mathfrak f} (t+\Delta t)\) before finding the state at \(\psi(t+\Delta t)\). To do this, we first find a predicted state \(\psi_a\) by assuming \(\Delta N_{\mathfrak f}=0\) and calculating \(\psi(t)\) according to the equation above. This lets us calculate an approximate \(\Delta N_{\mathfrak f} = N_{\mathfrak f a}  N_{\mathfrak f 0}\) and we use this in order to evolve \(\psi\). This is the ETD2RK scheme.
Note that \(N_{\mathfrak f}\) is a nonlinear function of the field variable \(\psi\), but can also be an explicit variable of time \(t\), i.e. \(N_{\mathfrak f}(\psi,t)\).
Therefore, in the code, it has to be encoded as a function of these two variables calc_nonlinear_evolution_function_f(self, psi, t)
.
For numerical purposes, it is useful to calculate the small \(\omega_{\mathfrak f}\) limit. We expand the exponential in its Taylor series and keep the leading order term to get:
In \(I_{\mathfrak f 1}\), and \(I_{\mathfrak f 2}\) there is a division by \(0\) when \(\omega_{\mathfrak f} = 0\). To avoid numerical issues related to this we use the above limits when \(\omega_{\mathfrak f}\) is smaller than a tolerance. We don't use the limit for \(I_{\mathfrak f 0}\) since it doesn't contain a division by \(0\). The function evolve_ETD2RK_loop
defined in the base system class performs an ETD2RK step. This function is called by the evolvers discussed in the model chapter if the method is defined as method = "ETD2RK"
. This is the default solver if method
is not set. The integrating factors for a given \(\omega_{\mathfrak f}(\mathbf{k})\) can be found with the function calc_evolution_integrating_factors_ETD2RK
where the variable tol
gives when the factors should be replaced by their leading order Taylor expansion.
Note that all solvers defined in the class \lstinline{BaseSystem} updates the time variable
self.t
to allow for timedependents in the nonlinear term.
The ETD4RK scheme
Following Ref.^{1}, we may generalize the method to a fourth order RungeKutta as follows
where
Algorithm: The ETD4RK scheme
In the small \(\omega_{\mathfrak f}\) limit, we have
Similar as for the EDT2RK case \(I_{\mathfrak f 1}\), \(I_{\mathfrak f 3}\), \(I_{\mathfrak f 4}\), and \(I_{\mathfrak f 5}\) contains a division by \(0\) when \(\omega_{\mathfrak f} = 0\).
We therfore replace these coeficients with their limits when \(\omega_{\mathfrak f}\) is smaller than a tolerance.
This has been important in order to make the the code stable for some of the systems.
In the same way as the EDT2RK scheme this is implemented as the function
self.evolve_ETD4RK_loop(self, integrating_factors_f, non_linear_evolution_function, field, field_f)
This function is called by the evolvers discussed in the model chapter if the method is defined as method = "ETD4RK"
, the integrating factors are found with
self.calc_evolution_integrating_factors_ETD4RK(self, omega_f, tol=10**(4))
.
The fully nonlinear limit
It is both interesting and enlightening to see the fully nonlinear limit of these equations, i.e., the limit in which \(\omega_{\mathfrak f} =0\), \(N_{\mathfrak f} = \partial_t \psi \equiv \dot{\psi}_{\mathfrak f}\) and the small \(\omega_{\mathfrak f}\) approximations are exact. For the ETD2RK scheme, we get
which is a twostage RungeKutta method called Heun's method.
The ETD4RK scheme becomes
Note that this is not the typical RungeKutta 4 method, due to the differences in calculating \(\psi_{\mathfrak f c}\). The reason is that a straightforward generalization of the RungeKutta 4 method will not produce a fourthorder method in the general case ^{1}.
The fully linear limit
If \(N=0\), the evolution equation changes to
An example of this is the schrödinger equation, for which \(\omega_{\mathfrak f} = \frac{1}{2} i \mathbf k^2\), so we get
This is an exact eqution, of course, so you may evolve this free particle solution to any time.
Testing
In order to test the numerical methods, study the simplest model of a field equation with a (nonlinear) forcing term, namely the heat equation
where \(T\) is the temperature in celcius, and \(f(\mathbf r)\) is a forcing term, which we model as
which represents a heating element with temperature \(T_0\) placed at \(\mathbf r_0\).
As a benchmark, we use the solve_ivp
of the scipy library sp.integrate
to solve the equation using a finite difference method.
The solutions match to a satisfactory degree, but a more thorough investigation into how the accuracy of the framework and integration methods scale with spatial and temporal resolution will be performed in the future.
Tests are included in test_base_system.py
, but for visual examination, here are animations of the initial condition \(T=0\) in all three dimensions
Algorithms for tracking defects
To be written
Calculating the velocity
The equations for the velocity are taken from Ref.^{2}, simplified using Mathematica and then substituted for python code using chatGPT.

Cox, S. M., & Matthews, P. C. (2002). Exponential Time Differencing for Stiff Systems. Journal of Computational Physics, 176(2), 430–455. https://doi.org/10.1006/jcph.2002.6995 ↩↩↩

Skogvoll, V., Rønning, J., Salvalaglio, M., & Angheluta, L. (2023). A unified field theory of topological defects and nonlinear local excitations. Npj Computational Materials, 9(1), Article 1. https://doi.org/10.1038/s41524023010776 ↩

Skogvoll, V., Angheluta, L., Skaugen, A., Salvalaglio, M., & Viñals, J. (2022). A phase field crystal theory of the kinematics of dislocation lines. Journal of the Mechanics and Physics of Solids, 166, 104932. https://doi.org/10.1016/j.jmps.2022.104932 ↩