The Vanadis model is a high-performance finite-element framework for three-dimensional atmospheric transport and dispersion. The spatial discretization is based on the Directional Residual (DR) method, which provides low-diffusion stabilization and good global mass-balance behavior in advection-dominated transport. Element matrices [H] and [C] are constructed on the CPU, while the solution stage employs a fully parallel Element-by-Element strategy without forming a conventional global sparse matrix. This structure maps efficiently onto GPU architectures: the CUDA implementation accelerates the solution of the linear system using a parallel Element-by-Element strategy combined with lightweight diagonal preconditioning.
Boundary conditions at the ground surface include either a Dirichlet constraint or a flux‑type deposition condition, with the latter controlled by the parameter α, which governs pollutant transfer to the surface.
Time integration is performed using an implicit time-discrete formulation, in which element matrices and right-hand-side contributions are constructed locally and the resulting algebraic system is solved iteratively at each time step. In the tested atmospheric dispersion cases, this approach demonstrated robust numerical stability and transient behavior for stiff advection–diffusion systems.
This solution is currently under patent evaluation. Patent application P.455475 (“pending”) has been filed with the Polish Patent Office. The full documentation will become publicly available once the Office completes the examination and publishes the application in its online database.
The source code of the Vanadis solver is distributed under a dual‑license model (GPL‑3.0 + Commercial License). The GitHub repository is archived and read‑only, serving as a permanent reference point. Users who wish to extend or experiment with the solver may do so by copying the ZIP archive into their own environment.
GIT repo: GitHub
Geometry: The model uses linear 8‑node hexahedral elements with standard shape functions Ni, which provide trilinear interpolation of geometry and field variables within each element.
Weighting functions: The weighting functions in the weak formulation are the Directional Residual (DR) functions, which provide a local, streamline-oriented stabilization of the advective transport operator. Unlike standard Galerkin weighting, the DR formulation modifies the test functions in the direction of advection while retaining the Element-by-Element structure of the discretization.
Shape functions (interpolation): The scalar field S(x,y,z) is interpolated using the standard trilinear basis
Linear equations:
Steady state solution:
Dynamic solution:
Time integration is performed using an implicit time-discrete formulation originating from numerical-methods lectures at AGH University of Science and Technology (Prof. Maciej Pietrzyk, ca. 1995). In the considered atmospheric dispersion simulations, the scheme demonstrated robust numerical stability and reliable iterative convergence for time steps of 100 s. Its ability to operate with relatively large time steps contributes substantially to the computational efficiency of large-scale transient advection-diffusion simulations.
.
This consistent transient-advection coupling is a fundamental part of the DR formulation rather than an implementation detail.
Vanadis runs on Linux or Windows systems and requires a Fortran compiler for building the CPU version of the solver. GPU acceleration is optional and uses the NVIDIA CUDA Toolkit.
For a representative test case consisting of approximately 2.8 million nodes and a 100-second time step, the following performance was observed:
CPU (Intel Core i5-4590 @ 3.30 GHz, 4 GB RAM): approximately 20 seconds to solve the linear system for a single time step.
GPU (NVIDIA GTX 1060, 3 GB): approximately 10 seconds to solve the linear system for a single time step.
These results illustrate the computational efficiency of the GPU-accelerated Element-by-Element solver compared with the single-core CPU implementation.
For the same large-grid case (approximately 2.8 million nodes), the Vanadis solver completes a 2,000-second transient simulation in approximately 20 minutes.
The results are visualized on two orthogonal surfaces: the ground surface and a vertical cross‑section. At time t=0 the first chimney is activated. After 1000 seconds the second chimney is switched on and the wind direction is changed. The influence of the tall building on the flow field and pollutant dispersion is clearly visible.
3D dynamic convection–advection model
Time step: 100 s
Pollution source: one emitter located at (x = 0 m, y = 0 m, h = 60 m)
Domain size: 8000 m × 8000 m × 420 m (height)
Example output: ground‑level concentration after 2000 s (shown on the left)
Data availability: full output for every time step on all visualization surfaces
Adjustable parameters:
Qv [kg/s] — emission rate
Kx, Ky, Kz [m²/s] — turbulent diffusion coefficients
P [1/s] — decay constant
vh [m/s] — vertical wind component
vx, vy [m/s] — horizontal wind components
α [m/s] — deposition / settling parameter
author: Marek Chodorski
e-mail: marek_ac@wp.pl