Numerical Simulations Laboratory

---

---

The featured projects presented below reflect our experience in scientific software development and numerical simulations for industry and academia. As software creators we are happy when our programs are used by many people and help in understanding complex phenomena in various fields of research.

Signal detection code for processing mass-spectrometry data

• Client: Mass Consortium Corporation, USA
• Technology: C++, OpenMP, Matlab
• Man hours: 280
• Scope: The developed C++ code for detection of features in mass-spectrometry data is based on CentWave algorithm described in a paper by Ralf Tautenhahn et al (BMC Bioinformatics 2008, 9:504, doi:10.1186/1471-2105-9-504). The locations and widths of chromatographic peaks are calculated entirely in Continious-Wavelet-Transform space by a simple and transparent procedure. The reading mass-spectrometry data files is done with help of two external libraries: MSToolkit and netCDF. Also, the LAPACK library is used in non-linear fitting of peaks data by Levenberg-Marquardt method. The code provides a fine control of the underlying algorithms by means of settings file, and it outputs detailed information about all stages of the detection procedure into special .features file. The output files and all data can be easily parsed and plotted by special Matlab scripts.
• Results: Due to improvements of the original algorithms the developed code performs better than famous xcms R package: it detects more peaks, and the peaks boundaries (in retention time) are calculated accurately. Typical results of chromatographic peaks detection in non-trivial case of overlapping peaks are demonstrated in the figure below. The left panels show the mass-spectrometry data in (m/z, time)-plane (top) and (intensity, time)-plane (bottom). Vertical red and blue lines mark maximums and boundaries (in retention time) of the detected peaks respectively. The right panels contain additional technical information about intermidiate steps of the algorithm.

Implementation of the boundary control method for inverse problem of acoustics

• Client: Saint Petersburg State University, Russia
• Technology: C++, OpenMP, MPI, Matlab
• Man hours: 800
• Scope: The Boundary Control Method (BCM) is ab initio approach to multidimensional inverse problems based on ideas and results of control theory, asymptotic methods (for propagation of discontinuities), functional analysis, and geometry. Its simple and clear mathematical background makes this rigourous method of rather general scope. In particular, the BCM is developed for acoustics, electrodynamics, Shrödinger equation, 1D Dirac system, 1D two-velocity system, reconstruction of Riemannian manifolds, metric graphs, and other dynamic systems. In this project we have implemented and thoroughly tested BCM for the sound speed determination in the acoustic equation on semiplane. The BCM provides time optimal step-by-step reconstruction procedure, which requires no ad hoc assumptions about the sound speed profile, and it works in the case of data given on a part of the domain boundary.
• Results: The application of the developed code to the realistic Marmousi model of acoustic medium demonstrates that the BCM is workable even in the case of extremely complicated and irregular field of acoustic rays. The results obtained open a broad perspective for using BCM in acoustic tomography. Find below the recently published papers about the BCM results (+ see references therein).

Fast universal hyperbolic solver (FUHS)

• Client: Saint Petersburg State University, Russia
• Technology: C++, OpenMP, MPI, Matlab
• Man hours: 1200
• Scope: The idea of the project is to develop an universal, efficient and easy-to-use solver for arbitrary system of hyperbolic equations in conservative form, so that any user with basic knowledge of C/C++ can implement a new set of equations within a working day. Such ambitious task can be achieved if central semi-discrete upwind third order numerical schemes of A. Kurganov are utilized. These finite volume truly multi-dimensional schemes: are developed for regular and irregular meshes, use weighted essentially non-oscillatory reconstruction, do not violate conservation laws, have low numerical diffision and dispersion. The code design is based on complete separation of mesh procedures, numerical scheme, system of PDEs, initial and boundary conditions, and the source terms. All these building blocks can be easily modified or replaced for any space dimension 1D, 2D or 3D. The job is distributed between computational nodes by means of MPI while OpenMP is used to parallelize calculations on each node. Simple input/output data formats and a set of ready-to-use Matlab scripts for data visualization make it possible to get first simulation results within a working day instead of months of coding and debugging.
• Results: The FUHS code is successfully used in numerous research projects.

Numerical study of flapping oscillations

• Client: Forschungszentrum Graz, Austria
• Technology: C++, OpenMP, MPI, Matlab
• Man hours: 200
• Scope: Kink-like magnetotail flapping oscillations in a Harris-like current sheet with earthward growing normal magnetic field component Bz are studied by means of time-dependent 2D linearized MHD numerical simulations. We use FUHS code with ∇ ⋅ B = 0 constraint enforced on each time step using the method of projection. For the first time a wave oscillating mode of the double-gradient magnetic instability has been observed and investigated. The dispersion relation and two-dimensional eigenfunctions are obtained. Coupled with previous results, present simulations confirm that the earthward/tailward growth direction of the Bz component acts as a switch between stable/unstable regimes of the flapping mode, while the mode dispersion curve is the same in both cases. It is confirmed that flapping oscillations may be triggered by a simple Gaussian initial perturbation of the Vz velocity.
• Results: Find attached the recently published paper (+ see references therein).

Shallow layer model of a dense gas dispersion

• Client: Software development company, Switzerland
• Technology: C++, Matlab
• Man hours: 400
• Scope: The purpose of the project is to develop a C++ code for simulations of a dense gas dispersion over complex terrains. The code is based on shallow layer approach developed in series of papers by Hankin et al [J. Hazard. Mater, A66, p. 211–261, 1999]. The model takes into account the effects of a gas cloud convection and buoyancy, the effect of ground slopes, surface and turbulent shear stresses, top entrainment and the leading edge terms that account for interaction among dense and ambient fluid. The shallow layer equations describing a gas cloud parameters (height, density, velocities) are solved by the flux-corrected transport algorithm of Zalesak, using corner transport scheme as a low order scheme and Zalesak fourth order scheme as a high order scheme. The code is implemented as a C++ library to be incorporated into the existing risk assessment software.
• Results: A typical simulation of a dense gas cloud propagation in a mountain region is presented below.

3D wind field model over complex terrains

• Client: Software development company, Switzerland
• Technology: C++, Matlab
• Man hours: 120
• Scope: The purpose of the project is to develop a C++ code for reconstruction of a 3D wind field from meteorological data in a given region with complex topography (mountains). We applied a modified version of a diagnostic wind model suggested in the paper by Ishikawa [J. of Applied Meteorology, V. 33, p. 733–743, 1994]. The approach ensures a divergence free wind field that has zero normal component at the ground surface and best fits the measured or the given wind parameters. Since a finite-difference discretization of equations leads to rather cumbersome matrix elements we have used Mathematica® to generate some parts of C++ sources (mathematical expressions) automatically to ensure that no errors are present in the code. Huge sparse linear system (matrix size ~ 200000 x 200000) is solved either by open-source (SuperLU) or commercial (Pardiso) library depending on compilation settings. The code is implemented as a C++ library to be incorporated into the existing risk assessment software.
• Results: Streamlines of typical reconstructed and adjusted wind fields are presented below.

Kinetic model of the interaction of magnetic perturbations with tokamak plasmas

• Client: Graz University of Technology, Austria
• Technology: C/C++, Fortran, Matlab, Mathematica
• Man hours: 2400
• Scope: The purpose of the project is to develop a kinetic linear model of the interaction of helical magnetic perturbations with cylindrical plasmas. The model is then implemented in the numerical code named KiLCA (Kinetic Linear Cylindrical Approximation). For a given equilibrium plasma configuration, KiLCA computes nonlocal plasma conductivity operator and solves Maxwell equations to calculate magnetic perturbation amplitudes inside the plasma. Derived quantities like the dissipated power density and the torques acting on the plasma are evaluated in postprocessing routines within KiLCA. A number of new and important results have been obtained by application of KiLCA to TEXTOR, JET and ITER tokamaks. Those results are well known by the fusion plasma community and are published in top level scientific journals.
• Results: Find attached the recently published papers (+ see references therein).

Upgrade of the KiLCA code to account for magnetic drifts

• Client: Graz University of Technology, Austria
• Technology: C/C++, Fortran, Matlab, Mathematica
• Man hours: 320
• Scope: The purpose of the project is to upgrade the existing KiLCA code taking into consideration (i) magnetic drifts of particles and (ii) energy conservation in particle collisions. The account of the magnetic drifts in nonlocal plasma conductivity operator is important near the so called "double null" point where a null of background electric field gets close to a resonant magnetic surface.
• Results: The improved plasma response model is used in self-consistent balance simulations to investigate the impact of external helical magnetic perturbations on evolution of plasma parameters in tokamaks.

Plasma formation during the laser ablation

• Client: Private person, Italy
• Technology: Fortran, Matlab
• Man hours: 160
• Scope: Laser ablation is used for a growing number of applications, such as pulsed laser deposition, nanoparticle manufacturing, micro-machining, surgery, chemical analysis and many others. Despite of active research for the last decade, the exact mechanisms are not yet fully understood due to many complex phenomena that occur and interact during the laser ablation process. Following a model proposed by A. Bogaerts et al. [Spectrochimica Acta Part B, 58, 1867–1893, 2003] we have developed a code that models the formation of a plasma near the target surface by numerical solving of the coupled Euler and Saha equations. The plasma code calculates the laser beam absorption and is coupled with an ANSYS solver that simulates the heat transfer in the target.
• Results: For high laser irradiances the plasma absorbs considerable amount of the beam energy and its formation has drastic impact on the laser ablation process.

Cross-equatorial abyssal currents in oceans

• Client: Private person, Canada
• Technology: Fortran, Matlab
• Man hours: 100
• Scope: Abyssal ocean currents, being denser than their surrounding fluid, are largely driven by gravity and steered by their local bathymetry. Over large scales they are also strongly influenced by the Coriolis force. The dynamics of these currents becomes complicated close to the equator, where geostrophic balance breaks down as the locally vertical component of the Earth’s rotation vector changes sign. Using a conventional approach based on the shallow water equations with the topography and the Coriolis force terms we have developed a code that models the cross-equatorial abyssal currents.
• Results: Here is a typical simulation of a cross-equatorial abyssal current. Near the equator y = 0 the Coriolis force cannot balance the force caused by the bottom topography slope in a parabolic channel (height(x,y) ~ x^2) and the current crosses the channel axis x = 0.

Ideal magnetohydrodynamics code for modeling the double-gradient instability

• Client: Forschungszentrum Graz, Austria
• Technology: C++, Matlab
• Man hours: 100
• Scope: The purpose of the project is to develop a code for detailed numerical investigation of the so called double-gradient magnetic instability, which is believed to be responsible for the magnetotail flapping oscillations - the fast vertical oscillations of the Earth’s magnetotail plasma sheet with a quasyperiod about 100-200 s. The set of linearized compressible ideal MHD equations in conservative form is solved by means of Lax-Friedrichs finite difference method. The results of numerical simulations confirm predictions of a simplified analytical theory. In particular, the instability growth rate is found to be close to a peak value provided by analytical estimates.
• Results: Find attached one of the recently published papers (+ see references therein).

Kink-like mode of a double gradient instability in a compressible plasma current sheet

Inverse code for reconstruction of time-varying reconnection rate and X-line location

• Client: Forschungszentrum Graz, Austria
• Technology: C++, Matlab
• Man hours: 220
• Scope: Magnetic reconnection is a physical process in highly conducting plasmas in which the magnetic topology is rearranged and magnetic energy is converted to kinetic energy, thermal energy, and particle acceleration. We have developed a numerical code for a remote-sensing method of reconstructing the reconnection rate and the location of X-line from single-spacecraft observations. The method is based on the two-dimensional analytical model of time-dependent Petschek-type magnetic reconnection in compressible plasma with asymmetric magnetic field configuration. The reconstruction technique has been successfully applied to series of flux transfer events recorded by spacecrafts in the Earth's magnetosphere.
• Results: Find attached one of the published papers (+ see references therein).

Reconstruction of the reconnection rate from Cluster measurements: Method improvements

The trolling calculator

• Client: Private person, USA
• Technology: Matlab
• Man hours: 24
• Scope: The purpose of this very small but also very interesting project is to predict the depth of an object towed through water at the end of a cable. The shape of the cable depends on boat velocity, properties of the object and the cable, temperature and salinity of water and is defined by the balance of tension, gravity, buoyancy, drag and lift forces. The resulting ODE is solved numerically and it allows to determine all parameters that are important for the successful trolling.
• Results: Here is the trolling calculator GUI.

---