Workshop on scientific computing 2016

June 9 - 12, 2016

Department of Mathematics, FNSPE CTU in Prague
Děčín, Czech Republic

List of abstracts

Spatially Isotropic Approximation Scheme for Motion by Mean Curvature
Michal Beneš, Jaromír Kukal
FNSPE CTU in Prague


The contribution summarizes the knowledge on the planar mean-curvature flow approximated by the Allen-Cahn equation with the gradient forcing  term. An overview on several spatial discretizations is presented including their advantages and drawbacks. Then a new approach respecting mesh directions but suppressing their influence on spatial isotropy of the solution is presented. Besides main idea, namely the obtained numerical results are presented and discussed.

Multidimensional exact solution for two-phase flow in porous media
Radek Fučík
FNSPE CTU in Prague


In general, analytical solutions serve a useful purpose to obtain better insights and to verify numerical codes. For flow of two incompressible and immiscible phases in homogeneous porous media without gravity, one such method that neglects capillary pressure in the solution was first developed by Buckley and Leverett (1942). Subsequently, McWhorter and Sunada (1990) derived an exact solution for the one and two dimensional cases that factored in capillary effects. This solution used a similarity transform that allowed to reduce the governing equations into a single ordinary differential equation (ODE) that can be further integrated into an equivalent integral equation. We present a revision to McWhorter and Sunada solution by extending the self-similar solution into a general multidimensional space. Inspired by the derivation proposed by McWhorter and Sunada (1990), we integrate the resulting ODE in the third and higher dimensions into a new integral equation that can be subsequently solved iteratively by means of numerical integration. We developed implementations of the iterative schemes for one- and higher dimensional cases that can be accessed online on the authors’ website.

A numerical method for optical flow
Viera Kleinová, Peter Frolkovič
SvF STU, Bratislava

Abstract: Optical flow is today very important topic in medicine, computer vision and image processing. The main goal is to determine optical flow based on level-set motion between two images. We present preliminary results of our numerical method. Our illustrative examples include synthetic and real data. Some representative results will be presented.

Multi-Material Remap for Staggered Arbitrary Lagrangian-Eulerian (ALE) Methods
Milan Kuchařík, Mikhail Shashkov
FNSPE CTU in Prague, XCP-4, LANL

Abstract: For hydrodynamic simulations of problems containing strong fluid compressions or expansions (such as simulations of laser-plasma interactions), Lagrangian methods employing a moving computational mesh are usually used. The Arbitrary Lagrangian-Eulerian (ALE) method is a successful approach preventing the mesh cells from numerical degeneracies. Such method typically consists of three steps: a Lagrangian solver, advancing the solution and the mesh in time; a mesh rezoner, keeping the moving mesh smooth; and a remapper, interpolating the fluid quantities between the meshes. Here, we focus on the last part of the ALE algorithm, especially on remapping of multi-material fluid quantities in the staggered discretization. In our remapping approach, all (both cell-centered and nodal) fluid quantities are remapped in a flux-form, while paying a special attention to their conservation and bound-preservation. Properties of our remapping method are demonstrated on a suite of selected hydrodynamic multi-material examples.

Useful Discrete Approximations of Laplacian and Gradient
Jaromír Kukal, Michal Beneš
FNSPE CTU in Prague, FNSPE CTU in Prague


Fourier and Laplace transforms were employed to develop novel difference scheme which is useful for the solving of parabolic partial differential equations and their systems. The first aim of scheme design is to develop such discrete approximation of Laplacian which has fixed accuracy order but has maximum possible radial order. When the grid is periodic, the coefficients of adequate scheme can be pre-calculated using linear methods. Tables of Laplacian coefficients are included for square, hexagonal, cubic and dodechedral grids. The second aim of scheme design is to develop adequate discrete approximation of gradient operator which indirectly satisfies the radial condition as follows. The gradiet operator has to generate Laplacian with maximum possible radial order. Unfortunately, the unknown coefficients of discrete gradient formula are roots of quadratic equation system. Finally, the tables of gradient coefficients are included for several particular cases together with the first simulation results.

TNL: FDM on GPU in C++
Tomáš Oberhuber
FNSPE CTU in Prague


We present Template Numerical Library with native support of CUDA for computations on GPUs. The library is written in C++ and it uses C++ templates extensively. The templated design of TNL allows to develop solvers of PDEs with GPU support relatively easily and almost without any knowledge of GPUs. The aim of TNL is to provide an easy to use tool for numerical mathematicians so that they may concentrate only to numerical methods but they can still profit from modern accelerators and parallel architectures. We will also discuss disadvantages of C++ templates and metaprogramming.

Numerical study of spiral motion and tip meandering
Petr Pauš, Shigetoshi Yazaki
FNSPE CTU in Prague, Meiji University, Tokyo


The talk focuses on the numerical simulation of spiral motion which occurs for example during Belousov-Zhabotinsky reaction. The spiral is simulated as an open parametric curve and approximated by the polygonal chain. The time evolution is based on the mean curvature flow equation. Parametric approach allows for the detailed description of the force applied to the spiral tip. The force consists of normal and tangential components which depend on several parameters and can be changed independently. We performed simulations under various settings of the tip force and studied the tip motion (meandering) in detail. The spiral motion is restricted to the circular domain. To avoid the loss of accuracy and expansion of the spiral outside the domain, we incorporated an algorithm which relocates the discretization points back to the computational domain allowing long time computations with relatively small number of discretization points.

Quantification of trapped gas in dual-porosity media with continuous and discontinuous domains
Michal Sněhota, Jan Šácha, Jan Hovind
FCE CTU in Prague, University Centre for Energy Efficient Buildings, CTU in Prague and Paul Scherrer Institut, Villigen, Switzerland


Nonwetting phase (residual air) is trapped in the porous media at water contents close to the saturation. Trapped gas phase resides in pores in form of bubbles, blobs or clusters forming residual gas saturation. In homogeneous soil media trapped gas is relatively stable until it is released upon porous media drainage. If porous media remain saturated, trapped gas can slowly dissolve in response to changed air solubility of surrounding water. In heterogeneous media, relatively rapid change in the trapped gas distribution can be observed soon after the gas is initially trapped during infiltration. It has been recently shown that the mass transfer of gas is directed from regions of fine porosity to regions of coarse porosity. The mass transfer was quantified by means of neutron tomography for the case of dual porosity sample under steady state flow. However the underlying mechanism of the gas mass transfer is still not clear. Based on the robust experience of visualization of the flow within heterogeneous samples, it seems that due to the huge local (microscopic) pressure gradients between contrasting pore radii the portion of faster flowing water becomes attracted into small pores of high capillary pressure. The process depends on the initial distribution of entrapped air which has to be considered as random in dependence on the history and circumstances of wetting/drying. In this study, the redistribution of trapped gas was quantitatively studied by 3D neutron imaging on samples composed of fine porous ceramic and coarse sand. The redistribution of water was studied under no-flow and steady state flow conditions. Two different inner geometries of the samples were developed. In the first case the low permeability regions (ceramics) were disconnected, while in the second structure, the fine porosity material was continuous from the top to the bottom of the sample. Quantitative 3D neutron tomography imaging revealed similar redistribution process in both cases of interconnected and disconnected fine pore systems. The rate of the redistribution was significantly higher in the case of steady state flow condition in comparison to no-flow conditions. The transfer from fine to large pores led to reduced hydraulic conductivity of the sample.

A Hybrid Parallel Numerical Algorithm for Three-Dimensional Phase Field Modeling of Crystal Growth
Pavel Strachota
FNSPE CTU in Prague


We introduce a hybrid parallel numerical algorithm for solving the phase field formulation of the anisotropic crystal growth during solidification. The implementation is based on the MPI and OpenMP standards. The algorithm has undergone a number of efficiency measurements and parallel profiling scenarios. We compare the results for several variants of the algorithm and decide on the most efficient solution.

Lattice Boltzmann Method and Natural Convection
Robert Straka
AGH Krakow

Abstract: The inclusion of a buoyant force into the LBM will be presented. Resulting macroscopic system of equations describe the problem of natural convection i.e. when the flow is induced due to density (temperature) gradients of given fluid. A good example is hot air movement above a heater. Multiple relaxation time (MRT) flavor of LBM together with Smagorinsky Subgrid Scale (SGS) LES turbulence model is applied to solve the fluid motion, single relaxation time (SRT) LBM for the second population of distribution functions, again with SGS LES is used to solve an advection-diffusion equation for the temperature field. Application of the above model is then applied to simulate heating of a room during the winter season. Dynamics of hot & cold air for different locations of the heater will be presented.

Two-stage Alpha-Stable Distribution Parameters Estimation Approach
Tran Quang Van
FNSPE CTU in Prague

Abstract: Financial asset returns tend to have heavier tail distribution than normal distribution and alpha stable distribution may be a suitable candidate for capturing this characteristic feature of asset returns. This heavy-tail distribution is characterized by four parameters which need to be estimated. They can be estimated by numerical integration approach, but it might be time consuming. We propose an approach based on maximum likelihood estimation method in which the parameters of alpha stable distribution are estimated consequentially in two stages. At the first stage two parameters alpha and beta are determined in an outer optimization loop from the standardized alpha stable distribution pdf obtained by fast Fourier transform. After that the remaining two parameters can be easily estimated in an inner optimization loop. The applicability of this two-phase likelihood maximization estimation technique is then verified on artificially simulated data alpha stably distributed and after that it is used to estimate parameters of alpha stable distribution of actual stock market indices returns series.

Construction of upper bounds of the HOMO-LUMO spectral gaps by semidefinite relaxation techniques
Daniel Ševčovič
Comenius University


In computational chemistry, the spectral properties of graphs describing organic molecules play an important role. The molecular orbital energy is associated with eigenvalues of the graph representing the molecule. More precisely, given an invertible graph G of an organic molecule, the energy of the highest occupied molecular orbital (HOMO) is related to the lowest positive eigenvalue of the adjacency matrix of the graph, and the energy of the lowest unoccupied molecular orbital (LUMO) corresponds to its largest negative eigenvalue. The so-called HOMO-LUMO spectral gap of the graph is then defined as the difference of HOMO and LUMO eigenvalues of the adjacency matrix of the structural graph of a molecule.

The aim of this talk is to present computational methods for obtaining upper estimates on the HOMO-LUMO spectral gap of graphs constructed from two prescribed structural graphs by bridging over their vertices. The problem leads to a mixed integer semidefinite programming problem which is an NP hard problem in general. We present a convex relaxation of the original problem leading to a numerically tractable  method for construction of the upper bound of the HOMO-LUMO spectral gap.

Efficient point cloud registration using PCL library
Róbert Špir, Lenka Hrapková, Karol Mikula
SvF STU, Bratislava

Abstract: Laser scanning and point cloud representation is common method to obtain 3D models of various objects or environments in medicine, architecture or digitalization of monuments and historical memorabilia. The most common problem is alignment of multiple scans of the same object from different angles to form single complete 3D model, especially when there is only small overlap of scans. In this work we will present efficient and fast registration of multiple large point cloud scans (more than 10 million points per scan) and their alignment using point-cloud library (PCL) with custom parallelization of calculation steps.