Student workshop on scientific computing 2016

June 9 - 12, 2016

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

List of abstracts

Mathematical modeling of phase interfaces of liquid mixtures using PC-SAFT equation of state.
David Celný
FNSPE CTU in Prague, Institute of Thermomechanics, Academy of Sciences of the Czech Republic


Phase interface is common phenomenon that can be found in each glass of mineral water. It is beneficial to understand its behavior and use it for more pressing matters such as carbon capture and storage technologies (CCS). Our work deals with gas-liquid phase interfaces of two component systems. These systems can form two basic types of phase interfaces on which we focus, namely planar interface (water level) and spherical interface (bubble). Both types of interface can be described via density function. One of the main tasks is obtaining molar density function profiles through minimization of grand potential functional. Our approach is based on theoretical results from Cahn-Hilliard gradient theory from 1957 that enables transformation of said problem into system of second order differential equation. Due to the performed simplification this system can be split into nonlinear algebraic system and one differential equation. Simplified formulation can be solved and provides resulting molar density function illustrating change of molar density over the change of coordinate (radius from center of bubble). Chosen approach requires the knowledge of mixture properties (chemical potential and pressure). These properties are computed with the PC-SAFT equation of state. We are able to compute the binary mixtures that exhibit polar interaction (i.e. CO2). Due to the unavailability of reference experimental data we investigated the system comprised of CO2and CH4.

On the correlation dimension estimation using rotational spectrum
Martin Dlask
FNSPE CTU in Prague

Abstract: Fractal sets are theoretical structures that don't exist in reality, however, many objects, such as biomedical or economic data have fractal pattern with noninteger dimension. Therefore it is necessary to develop methods that would provide unbiased estimation of fractal dimension. Technique of correlation sum belongs to the family of well-known and simple procedures which can be used for correlation dimension estimation, nevertheless, it was proven that the obtained dimension is underestimated. The presentation shows new methodology of correlation dimension estimation from point sets using rotational spectrum. Our novel approach employs the spectrum of point set which is averaged via rotation of power pattern. It was proven that the slope of the power spectrum can be stabilized by means of its rotation in a space with high dimension and the corresponding dimension estimate is unbiased. Precise correlation dimension estimation could be used later for data classification or prediction making.The efficiency of proposed method was tested using Monte Carlo simulation on the sets with known dimension.

Numerical analysis of the Lattice-Boltzmann method in 2D
Pavel Eichler
FNSPE CTU in Prague


In this contribution, numerical simulations based on the Lattice-Boltzmann method (LBM) will be discussed. First, the kinetic theory with the Boltzmann equation (BE) and some aspects of the derivation of the LBM from the BE with the BGK approximation of the collision operator will be introduced. In the last part, the application of LBM to the Hagen-Poiseuille flow between two infinite parallel plates will be discussed and the numerical results will be compared to the analytical solution in order to determine the experimental order of convergence.

Multicomponent transport model based on the Stefan-Maxwell theory and continuum mechanics
Jan Franců
FNSPE CTU in Prague

Abstract: Although the general equations had been proposed nearly 50 years ago, only recently the model combining the Stefan-Maxwell theory of diffusion with equations of fluid flow has been brought to light. In this contribution an elegant model for describing both convection and diffusion of components in a multicomponent fluid mixture is presented. Its effectiveness is demonstrated through a numerical analysis of simplified equations, which are solved using finite volume methods (FVM).

Template Numerical Library on Intel Xeon Phi
Vít Hanousek
FNSPE CTU in Prague


This talk will shortly introduce the Intel Xeon Phi coprocessor and its offload programing model. Then the current state of the Xeon Phi experimental support in the Template Numerical Library will be presented. Possible optimization techniques for a performance improvement will be discussed at the end of this talk.

Discontinuous Galerkin High-Order Nonlocal Transport and Energy Equations Scheme for Radiation-Hydrodynamics
Milan Holec
FNSPE CTU in Prague, Institute of Physics, AS CR, v.v.i

Abstract: The classical description of transport based on Chapman-Enskog approach has been always widely used in fluid models thanks to its simplicity. Nevertheless, it has been shown that the classical local approach is not accurate when the fluid parameters exhibit steep gradients, which is the typical case of laser heated plasmas. An intensive effort has been made to model the nonlocal radiative energy transport in radiation-hydrodynamics simulations in the last decades. From the existing models we solve directly the photon transport equation allowing one to take into account the effect of long-range photon transport. Our approach delivers a calculation efficiency and an inherent coupling of radiation to the fluid plasma parameters in an implicit way. The use of high-order discontinuous Galerkin method gives us an accurate solution to the transport, that obeys both limiting cases, i.e. the local diffusion asymptotic usually present in radiation hydrodynamics models and the collisionless transport asymptotic of free-streaming photons. In other words, we can analyze the radiation transport closure for the radiation-hydrodynamics and how it behaves when leaving the conditions of validity of Chapman-Enskog method. This is demonstrated numerically in the tests of the exact steady transport of any regime and the approximate time-dependent multi-group diffusion of energy. As an application we present simulation results of intense laser-target interaction, where the radiative energy transport, controlled by the mean free path of photons, shows the importance of the nonlocal model.

A massively parallel implementation of two-phase immiscible flow in porous media using the mixed-hybrid finite element method
Jakub Klinkovský
FNSPE CTU in Prague


The work deals with a numerical solution of two-phase immiscible flow in porous media and a massively parallel implementation of the model using the architecture of modern GPUs. We devise a semi-implicit numerical scheme that is based on the mixed-hybrid finite element and finite volume methods and stabilized using the upwind and mass-lumping methods. The scheme is implemented for parallel GPU architectures using the CUDA platform and the TNL library. The accuracy of the solver is verified by an experimental analysis of convergence

for benchmark problems with known semi-analytical solutions. For an advection-diffusion problem in heterogeneous porous medium, various capillarity models and numerical scheme variations are compared with a reference solution published in literature. The efficiency of the parallel computation on GPU is analyzed in detail for a selected test problem.

Numerical Solution of Constrained Curvature Driven Flow
Miroslav Kolář
FNSPE CTU in Prague


We investigate the numerical solution of the evolution law for the constrained curvature flow for open and closed curves in the plane. The model schematically reads as

normal velocity = curvature + force,

where the particular choice of the (possibly non-local) force term causes the structure-preserving property. In this contribution, we study the area preserving curvature flow, which originates in the theory of phase transitions for crystalline materials and originally describes the evolution of closed embedded curves with constant enclosed area. However, it can be also shown that this area preserving mechanism works for open curves with fixed endpoints as well.

The resulting motion law is treated by means of the parametric method, resulting in a system of degenerate parabolic partial differential equations. Unlike other possible approaches as, e.g., the level set method or the phase field method, the advantage of the direct approach is in the time efficiency and the ability to track the motion of open curves. 

We solve the parametric equations numerically by means of the semi-implicit flowing finite volume method. To enhance the numerical stability, we discuss the technique of tangential redistribution. Several results of our qualitative and quantitative computational studies will be presented.

Applications of planar and space curve evolution
Jiří Minarčík
FNSPE CTU in Prague

Abstract: In this contribution, we examine the theory of evolving curves and explore their use in several applications. We present a mathematical framework for describing curves in space which is a combination of the parametric and implicit approach. The framework has been developed to simulate the geodesic flow on stationary and moving surfaces. Both analytical and numerical results of the method will be presented. Along with applications in image processing, we will discuss the use of curves in modeling of the river channel centerline migration caused by the meandering process.

Mean Field Lévy Flight as Integer Optimization Heuristics
Matej Mojzeš
FNSPE CTU in Prague

Abstract: Integer optimization heuristics are the only feasible option for a variety of NP-hard optimization problems that need to be solved in real-world conditions. Based on trial and error and often enhanced by e.g. evolutionary, physical or biological processes they are able to find or approximate global optimum on very large search spaces. The purpose of our research is to contribute to family of heuristic method with a novel population based integer optimization heuristic that yields from the theory of Mean Field Annealing. Population center and covariance matrix are estimated for a given annealing temperature and then used as directional correction of Lévy Flight mutation, which delivers reputable results in heuristic optimization. Inspired by Competitive Differential Evolution, the proposed heuristic is of competitive nature with nine Lévy Flight mutations competing together and being selected according to their success. The resulting heuristic has four parameters: population size, regularization factor, annealing temperature and Lévy Flight temperature. Depending on the task complexity, there is relationship between searching efficiency and regularization, annealing and heavy-tailed flights. Last, but not least, performance of the novel method is demonstrated on benchmark Clerc's Zebra3 and Hilbert matrix inversion problems which are difficult tasks with many local extremes.

Numerical Computation of Two-Phase Compositional Flow in Porous Media
Ondřej Polívka
FNSPE CTU in Prague


We deal with the numerical modeling of compressible multicomponent two-phase flow in porous media with species transfer between the phases. 

The mathematical model is formulated by means of the extended Darcy's laws for all phases, components continuity equations, constitutive relations, and appropriate initial and boundary conditions. The splitting of components among the phases is described using a formulation of the local thermodynamic equilibrium which uses volume, temperature, and moles as specification variables.

The problem is solved numerically using a combination of the mixed-hybrid finite element method for the total flux discretization and the finite volume method for the discretization of continuity equations. These methods ensure the local mass balance. The resulting system of nonlinear algebraic equations is solved by the Newton-Raphson iterative method. The numerical flux is discretized in a way that no phase identification nor determination of correspondence between the phases on adjacent elements is required in contrast to the traditional approaches. This is very important for the simulations of CO2 sequestration because, typically, the CO2 is injected into a reservoir in the supercritical state at which the phase distinction is ambiguous. Moreover, our model performs well in situations where a phase appears or disappears and no switching of variables is needed.

We briefly describe the numerical method and provide several 2D simulations, e.g. CO2 injection into water saturated reservoir.

Computation of equilibrium states at constant internal energy, volume and moles
Tomáš Smejkal
FNSPE CTU in Prague

Abstract: In this contribution, phase stability and phase equilibrium of multicomponent mixtures at given internal energy, volume and moles will be discussed. We derive criterion for phase stability and devise numerical algorithm based on Newton-Raphson method for testing phase stability. We also devise a new algorithm for general p-phase equilibrium calculation, which is based on the direct maximization of the total entropy of the mixture with respect to the internal energy-, volume- and mole-balance constraints. We present the properties of the algorithms on several examples of phase equilibrium calculations.

Mathematical modeling of contaminant transport in porous media
Jakub Solovský
FNSPE CTU in Prague


This work deals with two phase compositional flow. We present equations describing two phase flow, component transport and interphase mass transfer. For this type of problems, we propose a numerical method based on the mixed hybrid finite element method. We implement several variations of this method using different approaches to solving resulting system of linear algebraic equations. We use direct and iterative solvers and parallel implementation using MPI. The method is verified on problems for which exact solutions are known or solutions can be found in literature. Numerical experiments show that the errors are similar for all variations of the method. The method is convergent and the order of convergence is slightly less than one. There are significant differences in the computational time. Iterative solvers are faster and the parallelism is advantageous while using fine meshes. In the next part, we focus on the compositional flow. Data from two experiments are used and numerical results are compared with measured values.

The first experiment was focused on evaporation of dissolved TCE. For low air flow rates above the water table, there is a good match with experimental data. For higher flow rates, the results differs. The second, more complex experiment in larger scale examined the influence of water table fluctuation and rainfall events on evaporation and transport of the dissolved TCE. For water table fluctuation, there is a good match with experimental data but for rainfall events there are significant differences. During the rainfall events there are uncertainties concerning the experimental data. Differences between the measured values and the numerical simulations indicate certain limits of the mathematical model used or the influence of other processes that are neglected in the current model. Finally we focus on hypothetical scenarios of vapor intrusion.

In the field scale, we examine effect of water table drop or rainfall events, that were in smaller scale studied experimentally.

Numerical results are similar to the second experiment.

A posteriori error estimates for finite element solutions of Poisson equation
Vojtěch Straka
FNSPE, CTU in Prague


A wide variety of phenomena in physics and other sciences can be described by partial differential equations. In majority of cases, finding an analytical (exact) solution is not possible, therefore numerical methods are applied as simulation tools. However, these methods typically only deliver an approximate solution, which is different from the exact solution. For evaluation of the error between the known numerical and the unknown exact solution, a posteriori error estimates can be used. 

In this presentation, a general introduction to a posteriori error estimation will be made. Then a specific form of a posteriori estimates for Poisson equation will be discussed and finally numerical results for model problems will be presented.

Quantitative evaluation of water distribution from two and three-dimensional neutron images the during ponded infiltration
Jan Šácha, Michal Sněhota, Jan Hovind
FCE CTU in Prague, Paul Scherrer Institute, Switzerland


Modern imaging techniques such as neutron imaging (NI) provide spatial and temporal information about the water and air distribution within the porous media. This information during hydrological processes is important for evaluating current and developing new water transport models. NI is characterized by relatively short acquisition times (seconds) and high resolution of images (micrometers). The acquisition time increases with the better resolution and vice versa. Depending on a research focus (static or dynamic processes) the choice of parameters is of a high importance. At the same time, the appropriate data processing has to be applied to obtain results free of bias and artifacts. Ponded infiltration experiments were conducted on two soil samples packed into the quartz glass columns of inner diameter of 29 and 34 mm, respectively. First sample was prepared by packing of fine and coarse fractions of sand and the second sample was packed using coarse sand and disks of fine porous ceramic. Ponded infiltration experiments conducted on both samples were monitored by neutron radiography to produce two dimensional (2D radiograms) projection of images during the transient phase of infiltration. During the steady state flow stage of experiments neutron tomography was utilized to obtain three-dimensional (3D tomograms) information on gradual water redistribution. The acquired radiographic images were normalized for background noise and spatial inhomogeneity of the detector, fluctuations of the neutron flux in time and for spatial inhomogeneity of the neutron beam. The radiograms of dry sample were subtracted from all subsequent radiograms to determine water thickness in the 2D projection images.

All projections were corrected for beam hardening and neutron scattering by empirical method. Parameters of the correction method uses were identified by fitting the volume of water in the entire sample in given time (from radiograms) to gravimetrically determined amount of water in the sample. The results from this correction is 2D water thickness maps of the sample. Tomography images were reconstructed from corrected water thickness maps to obtain the 3D spatial distribution of water content within the sample which can be compared with results of mathematical models.

Modeling of Water-Ice Interface within Freezing Soil at Micro-Scale
Alexandr Žák
FNSPE CTU in Prague


This contribution deals with a 2D micro-scale model of thermomechanical processes during solidification of a medium within porous material. The problem description is performed by means of heat balance and momentum conservation within individual phases; the solidification is traced using a phase-field equation. Suitable couplings of multi-phase and multi-physics are introduced. Several computational results are