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 advection and strict mass conservation, making it suitable as a reference solution for evaluating stabilized schemes. Element matrices [H] and [C] are assembled on the CPU, while the solution stage employs a fully parallel element‑by‑element strategy without forming a global 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 inside the solver, taking advantage of the diagonal structure preserved by DR stabilization.
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 penetration into the surface.
Time integration is performed using a custom stabilized implicit second‑order multistep scheme inspired by BDF‑type methods. This approach provides strong numerical stability and accurate transient behavior for stiff advection–diffusion systems. The combination of DR spatial discretization, GPU‑accelerated element‑wise solving, and the stabilized multistep time scheme results in a robust and efficient tool for large‑scale atmospheric dispersion simulations.
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.
DOI: 10.5281/zenodo.21352014
GIT repo: https://github.com/marek-ac/vanadis_solver-directional_residual/releases/tag/v2026.1.1
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, ensuring consistency between the residual minimization and the spatial discretization. Unlike Galerkin weighting, DR preserves the diagonal structure required for the element‑by‑element solver.
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 a stabilized implicit second-order multistep scheme originating from AGH numerical-methods lectures (prof. Maciej Pietrzyk, ca. 1995), which provides robust stability for stiff transport problems and ensures accurate temporal evolution even under rapidly changing flow conditions. This choice of integrator enables the solver to maintain high accuracy without requiring excessively small time steps.
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 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 to the single‑core CPU implementation.
For the large grid (≈ 2.8 M nodes), the DR formulation completes the full 2 000‑second transient simulation in 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