GitHub

Release notes

The format used for this changelog is based on Keep a Changelog, and this project adheres to Semantic Versioning. Notice that until the package reaches version v1.0.0 minor releases are likely to be breaking. Starting from version v0.3.1 breaking changes will be recorded here.

Unreleased

Added

  • AMG() is a preconditioner, SolverSetup(solver=Cg(), preconditioner=AMG()), with Cg(), Bicgstab() or the new Richardson(). AMG() maps to GAMG() on distributed meshes, and GAMG()/BoomerAMG() to the native AMG on serial ones. #244
  • Added examples/transient_verification/, six small cases comparing Euler with CrankNicolson against exact or reference solutions and printing the observed temporal order. #228
  • Runtime post-processing (postprocess) runs on distributed meshes in piso!, simple! and laplace!, and laplace! post-processed scalars are written to OpenFOAM output with zero-gradient patches. #234
  • Distributed gate tests share one mpiexec session, and rank counts above XCALIBRE_TEST_MAX_RANKS (default 2) are skipped. #225
  • Added the test grid test/grids/OF_graded_channel, a 700-cell 2-D OpenFOAM channel with streamwise grading for transient tests. #237
  • The density-based solver godunov! runs on distributed meshes with all flux, reconstruction and time-stepping schemes and the Laminar, KOmega and KOmegaSST models. #222
  • Added a characteristic far-field treatment in godunov! for patches where U, p and T are all Dirichlet or DirichletFunction, so acoustic waves leave the domain. #217
  • Added local pseudo-time stepping Runtime(adaptive=LocalTimeStepping(maxCo, startCo, ramp)) and the Jameson multistage schemes RK3(), RK4() and RK5(), with rhoU and rhoE residuals and force histories. #216
  • Added the LeastSquares gradient scheme, selected with Schemes(gradient=...), which is exact for linear fields on skewed meshes; it is off by default. #215
  • Added consistent (SIMPLEC) coupling with run!(model, config; consistent=true) for simple! and csimple!; it is off by default. #147
  • Added the second-order upwind divergence schemes LinearUpwind and LimitedLinearUpwind, which keep the matrix diagonally dominant; other equations fall back to Upwind with a warning. #169
  • he can be convected with LinearUpwind and LimitedLinearUpwind, instead of stopping with an ArgumentError. #207
  • Added the divergence scheme LUSTDC, a deferred-correction form of LUST for steady runs that do not converge with LUST; LUST is unchanged. #211
  • Added the divergence scheme LUSTLU, an opt-in second-order variant of LUST for transient and LES runs; it is not available for he. #213
  • Added the Venkatakrishnan gradient limiter, Schemes(limiter=Venkatakrishnan(K)), next to CellBased. #206
  • Added reorder_mesh!, which renumbers cells, faces and nodes for memory locality; it is opt-in and mesh files are not touched. #188
  • Added mixed-precision linear solves with SolverSetup(...; precision=MixedF32()) and MixedF16(), which solve corrections in lower precision while the fields stay in full precision. #189
  • activate_multithread(CPU(static=true)) and first_touch(mesh) place memory by parallel first touch on NUMA nodes; call mesh = first_touch(mesh) before building the model. #167, #190
  • Added volume_integral, weighted_volume_integral, volume_average and total_volume for backend-agnostic volume integration. #137
  • Added the Robin boundary condition for mixed constraints a·φ + b·∇φ·n = value. #138
  • Added the WALE LES model, LES{WALE}(), which also runs on distributed meshes. #177
  • run! takes a progress keyword (default true); set progress=false to hide the progress bar, e.g. for batch jobs. #162
  • Added distributed (MPI) simulation: distribute gives each rank its partition from the same serial case setup, with PETSc linear solves and OpenFOAM decomposed output. #162
  • Added partition_mesh for offline decomposition and distribute(dir) to load one part per rank. #162
  • Added restart for distributed incompressible runs with run!(...; restart=...), using OpenFOAM binary output. #162
  • Added distribute(FOAMCase(dir)) to read a decomposePar case, and repartition to rebalance a distributed mesh. #162
  • Added is_root, which is true on rank 0 and in serial runs. #162
  • Added BoomerAMG and GAMG pressure preconditioners for distributed solves, and petsc_options on run!. #162
  • Distributed solves map Cgs to PETSc cgs, and DILU, ILU0GPU and IC0GPU to per-rank block ILU(0) or ICC(0) with a warning; PETSc.set_petscint!(Int32) selects 32-bit indices. #162
  • Distributed GPU runs detect whether MPI is CUDA-aware and stage through host memory with a warning when it is not. #162
  • Added the wall_distance keyword for KOmegaSST and KOmegaLKE, with MeshWave() (the new default) or Poisson(). #170
  • Added div!(phi::VectorField, Γf, tensor, BCs, config), the divergence of a cell tensor scaled by a face coefficient, and the lazy Dev2(A) tensor. #168

Changed

  • Distributed Bicgstab() and Gmres() solves with GAMG(), BoomerAMG() or IC0GPU() no longer store a second, scaled copy of the matrix, lowering peak memory; the scaled matrix is applied on the fly. #250
  • With an AMG preconditioner, equations no longer build an unused second preconditioner; itmax defaults to 1000 for every solver. #244
  • The residual of each field is evaluated on its value before the linear solve, so residual histories no longer depend on the linear solver, its precision or tolerances, or on serial, distributed or GPU runs. Histories shift by one iteration, and a field that starts at zero reports the unnormalised Σ|b - Ax| on its first iteration. #243
  • atol and rtol mean the same on every linear-solve path: each solve stops when the Jacobi correction ‖D⁻¹(b − Ax)‖₂ ≤ atol + rtol·‖D⁻¹(b − Ax₀)‖₂, so atol is in the field's units whatever the solver, preconditioner, precision or backend, and warnings print that norm. Jacobi solves keep their tolerances; Cg() stops tighter than before (it used √(rᵀD⁻¹r)), Gmres() restarts every 30 iterations, and PETSc adds atol and rtol instead of taking their maximum. #245
  • The CI and documentation workflows cancel superseded pull request runs and queue runs for main and tags. #202
  • The explicit viscous stress and div! of cell tensors run faster on GPUs, with unchanged results. #187
  • xcal_foreach on CPU(static=true) uses fixed chunks matching the threaded SpMV and vector operations. #167
  • Threaded loops start their chunk tasks as a tree, which speeds up many-thread solves; results are unchanged. #191
  • activate_multithread(backend::CPU) defaults nthreads back to 1, as Krylov vector operations use Julia threads rather than BLAS. #162
  • The multiphase and thin-film solvers no longer print an @time summary of their iteration loop. #162
  • run! with a negative write_interval builds no output writer, reducing run start-up cost. #162
  • Kernels compile once for any mesh size, patch count or rank count. #162
  • The minimum supported Julia version is 1.10, and CI runs on Julia 1.10 and 1.13. #162, #165
  • Distributed VTK output stops with an error; OpenFOAM fields carry dimensions entries. #162
  • Mesh2 and Mesh3 store each property as its own array field; mesh.cells, mesh.faces and mesh.nodes still index as before. #162
  • Custom @define_boundary bodies and scheme! methods receive cells and faces arrays instead of one cell and face. #162
  • Discretisation kernels receive equation terms and arrays without the mesh, reducing GPU per-thread memory. #162
  • On the CPU, Krylov vector operations, SpMV and Jacobi run on Julia threads over a fixed row split. #162
  • Periodic boundary conditions keep an Int32 mesh's equation matrix at Int32. #162
  • The AMG hierarchy uses the index type of the equation matrix. #162
  • Global cell and face ids of distributed meshes are Int64. #162
  • The explicit viscous stress takes boundary-face gradients from the boundary condition. #181, #186

Fixed

  • The boundary face delta, and with it the wall-function distance wall_normal_distance, is the distance from the owner cell centre to the nearest point of the face, set once when the mesh geometry is built by FOAM3D_mesh, UNV3D_mesh, UNV2D_mesh and the 2D block mesher. It equals (Cf - C)·n when the normal from the cell centre meets the face, and on skewed near-wall cells it no longer collapses towards zero, which made y+ vanish and the wall value of omega grow by orders of magnitude. #247
  • Float32 GAMG() on distributed matrices above 2^24 + 1 global rows now stops at solver construction with an error naming FullPrecision(), instead of failing inside PETSc setup. #254
  • UNV2D_mesh cell centres are the area centroids of the cells, and the face delta, e and interpolation weights are computed from these final centres. Previously the face geometry was built from the vertex-average centres, which were then replaced by a perimeter-weighted average, so the stored face geometry did not match the cell centres on triangular and irregular cells. #251
  • Every dependency now has a [compat] entry, including the standard libraries LinearAlgebra, Logging, Printf and SparseArrays, as the General registry requires. #242
  • Documentation: the release notes link to the correct repository, and Robin, volume_integral, weighted_volume_integral, volume_average and total_volume are now documented. #242
  • ConstantViscosity is initialised correctly and its test is registered in the suite. #204
  • Distributed linear solves with PETSc.jl 0.5 no longer print its own convergence warning next to the XCALibre.jl one; unconverged solves warn once. #205
  • Tests: test/distributed/test_perf.jl now checks the correct halo-exchange budget and converged solves, and runs in the two-rank gate session #240
  • Distributed linear solves that stop without meeting atol or rtol now warn with the method, equation, reason and residual, at most 10 times per equation; converged solves print nothing. #239
  • Linear solves with Bicgstab() or Cgs() no longer return NaN on breakdown; a broken-down solve is repeated from its initial guess with a fixed shadow vector, and warns if it breaks down again. #238
  • Results no longer change from run to run with thread scheduling: boundary-face contributions are summed per owner cell in face order, so two identical runs on the same thread count are bit for bit identical. #232
  • bounding_box(mesh) is a reduction over the mesh nodes returning SVectors, and the Empty distance is the node extent along the patch normal, fixing crashes and wrong LES filter widths. #233
  • The thin-film solver (Momentum{EFM}, filmModel!) now gives a steady state that does not depend on the time step, with prescribed inflow fluxes honoured exactly. #231
  • Qcriterion now honours its start, stop and update_interval keywords, computing Q only inside its window. #229
  • The Laplace example and tests pass Schemes(...) and SolverSetup(...) directly instead of a named tuple. #235
  • The kinetic-energy face interpolation of the enthalpy equation works on KernelAbstractions 0.9.25 and later, as the compat range allows. #236
  • viscous_force and wall_shear_stress work on DirichletFunction, RotatingWall and zero-gradient patches, and viscous_force, wall_shear_stress and stress_tensor take a time keyword. #230
  • The explicit viscous stress and viscous work use the boundary-face velocity gradient of every velocity condition, and periodic faces are treated conservatively. #227
  • CrankNicolson is now the second-order trapezoidal rule inside Time, with CrankNicolson(ψ) setting the off-centring; it raises an ArgumentError in multiphase!, filmModel!, godunov! and steady solvers. Euler and steady results are unchanged. #226
  • The periodic face-flux correction uses the face coefficient of the periodic Laplacian, so continuity holds where the partner cell is also a neighbour; other periodic incompressible results are unchanged. #226
  • Transient piso!, cpiso! and multiphase! runs reach the same steady state for different time steps, with a bounded blend of the old face flux. Euler results change slightly; temporal order is unchanged. #225
  • godunov! with LocalTimeStepping compiles on GPUs, as turbulence_dt is now a singleton type; the :min, :mean and :max inputs are unchanged. #218
  • pressure_force and viscous_force return the whole-patch force on every rank of a distributed mesh. #219
  • Ghost cells of bounded turbulence fields hold the owners' values on distributed meshes, as they are exchanged after clipping. #220
  • output=OpenFOAM() writes fields without their own boundary conditions, such as rho and y, with zero-gradient patches. #221
  • viscous_force and the forces= history of godunov! return zero viscous force on Slip and Symmetry patches. #223
  • godunov! runs the KOmega, KOmegaSST and LKE RANS models, which previously failed at initialisation. #214
  • Documentation: the numerical setup page no longer shows the divergence-scheme text in the navigation bar, and describes the SupersonicFlow options. #209.
  • godunov! defaults reconstruction to Upwind() when it is not given. #210.
  • Fixed the Smagorinsky, WALE and KEquation LES models failing on GPUs with "Scalar indexing is disallowed". #212
  • Fixed the sensible-enthalpy equation creating total enthalpy when he uses a scheme other than Upwind or BoundedUpwind; schemes without a defined face value raise an ArgumentError. #203
  • Fixed the face mass flux of simple!, piso!, csimple!, cpiso!, the MRF simple! and potential_flow! not conserving mass when ncorrectors > 0; results of such runs change. #199
  • Fixed the viscous work of the sensible-enthalpy equation of csimple! and cpiso! creating energy in wall cells. #198
  • Fixed the examples/2D_cylinder_heated.jl example failing with UndefKeywordError, as cp is now supplied. #200
  • Fixed continued runs not retracing the uninterrupted run, as the last face mass flux is now reused. #196
  • Fixed the distributed continuation test failing under multithreading; residual histories are compared to 1e-5. #201
  • Fixed the face mass flux of csimple! and cpiso! for Compressible fluids differing from the flux the pressure equation balances. #195
  • Fixed the sensible-enthalpy equation of csimple! and cpiso! counting viscous heating twice. #194
  • Fixed compressible solvers evaluating energy boundary values at the wrong time. #192
  • Fixed FOAM3D_mesh failing on meshes whose neighbour list pads boundary faces with -1. #193
  • Fixed the explicit viscous stress destabilising turbulent runs on skewed near-wall cells; turbulent results change slightly. #180
  • Fixed the incompressible momentum equation missing the dev2 part of the viscous stress, now included as an explicit source. #168
  • Fixed warm-started serial Bicgstab solves degrading as steady runs converge. #166
  • Fixed bound! giving schedule-dependent results, now deterministic on every backend. #174
  • Fixed the KOmegaLKE diffusion terms, which now include the intermittency factor and use boundary values of the fields. #173
  • Fixed the KOmegaLKE production term, now the double contraction of the velocity gradient; results in flows with shear and strain change. #172
  • Fixed Neumann Laplacian{Linear} ignoring the sign of the Laplacian term. #164
  • Fixed CPU linear solves failing when called inside threaded loops or concurrent tasks. #162
  • Fixed AutoTune() launches over an empty range failing. #162
  • Fixed distributed runs hanging when one rank fails; the error is now raised on every rank. #162
  • Fixed distribute(reader; dir) silently reusing another mesh's parts; key repartitions when it differs. #162
  • Fixed mesh readers failing or writing out of bounds on connectivity overflow; the error asks for integer_type=Int64. #162
  • Fixed an intermittent segmentation fault in distributed GPU runs. #162
  • Fixed turbulence models computing in Float64 on Float32 meshes, which now use the mesh float type. #176
  • Fixed wall functions using the floored boundary distance; they now use wall_normal_distance. #183
  • Fixed the AMG zero-allocation coarse-solve test failing intermittently. #185
  • viscous_force and wall_shear_stress integrate the boundary viscous stress the momentum equation applies, -ρ ν_eff (∇U + ∇Uᵀ - ⅔ tr(∇U) I)⋅n with the boundary-face velocity gradient, instead of its tangential normal-gradient part only. The transposed and normal stress terms now enter the reported forces, which therefore change; zero-gradient patches report the force of their transposed stress, while Slip and Symmetry patches still report none. #249

Breaking

  • CrankNicolson is the trapezoidal rule inside Time and raises an ArgumentError in multiphase!, filmModel!, godunov! and steady solvers. #226
  • Results change by design with the corrected mass flux of runs with ncorrectors > 0 #199, transient piso!, cpiso! and multiphase! #225, compressible viscous heating #194, KOmegaLKE production #172, and the explicit viscous stress #180.
  • Turbulence models on Float32 meshes compute in Float32, so their results change. #176
  • Distributed VTK output raises an error. #162
  • activate_multithread(backend::CPU) defaults nthreads to 1. #162
  • potential_flow! starts from rest by default (from_rest=true); from_rest=false keeps the previous behaviour. #182
  • KOmegaSST and KOmegaLKE default to the MeshWave() wall distance, changing their results; wall_distance=Poisson() restores the previous method. #170
  • setField_Box!, setField_Circle2D!, setField_Sphere3D! and setField_Expression! take a required config keyword and run on CPU and GPU. #137
  • Hand-written boundary-condition functors and scheme! methods take cells and faces and must index them. #162
  • The Mesh2 and Mesh3 type parameters changed, so code dispatching on them must be updated; mesh.cell_nsign is Int8. #162
  • PETSc.jl 0.4.14 and 0.5 are both supported, resolved by Julia version. #165
  • MPI and Metis are hard dependencies. #162
  • Mesh readers default to integer_type=Int32; meshes exceeding 2^31 need integer_type=Int64. #162

Deprecated

  • SolverSetup(solver=AMG(mode=...)) runs as solver=<mode>, preconditioner=AMG(...); AMGSolver() is now Richardson(). #244

Removed

  • Removed the unused xmul and xmul! methods. #162

Version v0.6.1 - 2026-09-25

Added

  • Added FixedHeatFlux boundary condition #149
  • Added AMG-preconditioned stabilized biconjugate gradient for solving non-symmetric equations/ #150
  • Added potential_flow! to initialise a simulation from a divergence-free potential-flow field. Velocity boundary conditions supply the initial face flux, and velocity-potential boundary conditions are inferred from the pressure boundary conditions: fixed pressure becomes fixed zero potential, periodic patches stay periodic, and every other patch uses zero normal gradient. Supports non-orthogonal correctors through ncorrectors #158

Changed

  • 3D mesh geometry is now computed by a single shared routine used by both the UNV3 and OpenFOAM readers, which previously disagreed on cell centroids and volumes for the same mesh #153
  • The OpenFOAM writer now preserves a complete existing constant/polyMesh instead of overwriting it, so results can be written alongside the original mesh #153
  • The OpenFOAM writer compares the point and face counts declared by an existing constant/polyMesh with the simulation mesh, and warns and rewrites the mesh files when they disagree #153
  • 2D_compression_corner.jl and 2D_cylinder_transonic_RANS.jl now use Slip on patches that are physically symmetry planes or inviscid walls, replacing Zerogradient #156
  • _nonorthogonal_face_correction now builds the correction from a projected-delta correction vector rather than the over-relaxed Ef/T_hat decomposition, and correction_weight is no longer needed #158
  • CSIMPLE now relaxes pressure explicitly for Compressible fluids, as it already did for WeaklyCompressible ones, instead of relaxing the pressure equation implicitly through irelax. Relaxation had to move after the mass-flux correction so that the flux is built from the unrelaxed solution, and implicit relaxation cannot be deferred in the same way because it alters the matrix the flux is reconstructed from. Compressible cases converge along a different path and may need their pressure relaxation factor retuned #158
  • 2D_bump_NASA_kwSST.jl, 2D_compression_corner.jl and 2D_cylinder_transonic_RANS.jl solver tolerances and relaxation factors retuned for the corrected fluxes #158
  • Cached the sparse index maps the assembly writes through. The CSR sparsity pattern is fixed for the life of an equation, so ScalarEquation and VectorEquation now carry diag_nz and face_nz, the nzval index of the diagonal and of each cell-neighbour coefficient, resolved once at construction instead of searched for on every assembly. Both types gain a type parameter and two fields #161
  • Cached the Laplacian face coefficient on the mesh. Mesh2 and Mesh3 gain a type parameter and the face_gDiff field, derived from the face geometry at construction rather than passed in, and refreshed by update_face_gDiff! for readers that fill the geometry afterwards. The internal-face form norm(((Sf⋅Sf)/(Sf⋅e))*e)/delta reduces exactly to area/(|normal⋅e|*delta), and boundary faces keep area/delta, which is what every @define_boundary Laplacian block already uses. A mesh serialised with an earlier version can no longer be loaded and must be regenerated #161
  • Turbulence models now build the wall-function accumulators once in initialise instead of allocating and zeroing two cell-sized arrays on every outer iteration. The pair is nothing when no patch uses a wall function #161
  • Threaded the diagonal preconditioner application. LinearOperators' opDiagonal applies the diagonal with a serial broadcast, once (Cg) or twice (Bicgstab) per Krylov iteration, and threading it is worth about 6% of every Krylov iteration on the motorBike case #161
  • The SIMPLE and KOmega flux fields are built without a mesh reference, and their consumers take the mesh from the cell field. A field struct that stores a mesh costs 416 bytes of per-thread local memory (an adapted Mesh3) in every GPU kernel whose argument aggregate is not scalarised, which is 13.6% of a motorBike GPU run #161
  • The KOmega cell and face loops keep the mesh out of the kernel. xcal_foreach dispatches on the backend so the GPU path calls the closure directly, where AcceleratedKernels passed it to a callee that is not inlined and spilled the captured environment, and both loops bind their fields to their bare storage. The production loop drops from 7112 bytes of stack frame and 11.0 ms per launch to none and 0.20 ms, the viscosity loop from 1840 bytes and 8.86 ms to none and 0.12 ms. On the CPU the AcceleratedKernels task partitioning is unchanged. Together with the field change above, motorBike runs 500 iterations in 21.0 s against 26.3 s #161
  • Removed the unused kf and omegaf fields from KOmega. They were allocated and destructured but never read, and were carried into every kernel that takes the turbulence model #161
  • Removed the undocumented XCALIBRE_AMG_DEVICE_COARSE_SOLVE, XCALIBRE_AMG_DEVICE_COARSE_MAXITER and XCALIBRE_AMG_DEVICE_COARSE_RTOL environment variables. They enabled a second, parallel device coarse-solve path that duplicated coarse_solve = OnDeviceKrylov() and reached it only with mode = Cg(), the one combination the AMG constructor rejects because an inner Krylov solve is nonlinear in its right-hand side and breaks the fixed-operator assumption of preconditioned CG. The documented coarse_solve keyword is now the only way to choose how the coarsest level is solved, and it is described in the user guide #159
  • Removed the coarse_tmp scratch vector from AMGLevel. It was written only by the removed path, so every level allocated a full-length vector on the device that nothing read #159

Fixed

  • Fixed mesh readers orienting some face normals into their owner cell. The orientation was decided by testing each normal against an estimated cell centre, which lies on the wrong side of a face in skewed or concave cells: 39 cells of a 0.6M-cell and 92 of a 10M-cell snappyHexMesh motorBike mesh did not close (sum of outward face areas up to 52% of the cell surface), and a concave 2D cell got wrong nsign entries and a wrong area. Faces are now oriented from each cell's topology (edges traversed in opposite directions by adjacent faces, overall sign from the divergence theorem), which is exact for any closed cell, in UNV3D_mesh, UNV2D_mesh and BlockMesher2D; FOAM3D_mesh keeps OpenFOAM's own orientation, which the format defines. nsign is now set from ownership in 2D, as in 3D, and 2D cell areas no longer depend on the orientation test. The VTK writer orders polyhedral face nodes from the same orientation instead of its own centre test #163
  • Fixed missing density term in set_production! needed for cases with non-unity density #151
  • Fixed wall_shear_stress to apply the effective viscosity (nueff) scaling to the x-component of the shear stress vector - previously only the y and z components were scaled #152
  • Fixed OpenFOAM boundary file parsing so that patch groups (inGroups) no longer corrupt the patch list. The parser is now token based, ignores unknown dictionary entries, and reports malformed files with an ArgumentError #153
  • Fixed OpenFOAM mesh writer truncating point coordinates to six significant figures, which prevented a written mesh from round-tripping #153
  • Fixed reconstruct! omitting boundary faces when rebuilding a cell vector from a face-normal flux. The accumulation looped over cell_faces, which holds internal faces only. A cell becomes genuinely rank deficient when a direction is spanned by no internal face at all, as on a one cell thick mesh with an empty patch: on OF_pitzDaily the through-thickness moment is zero for all 12225 cells and the reconstruction returns the zero vector everywhere. Boundary faces now contribute to the least-squares system #155
  • Fixed the invertibility test in reconstruct! comparing an area-cubed determinant against the absolute tolerance eps(TF). The moments scale with face area, so on a millimetre-scale mesh the determinant of a perfectly conditioned interior cell is already below eps (6.4e-29 on the 1 mm lid-driven cavity) and the reconstruction is silently zeroed everywhere, interior cells included, and at far coarser scales in Float32. The test is now scaled by the magnitude of the moments #155
  • Fixed Wall and RotatingWall diffusion imposing only the tangential part of the boundary velocity together with the cell's own normal component, which let the wall-normal velocity float instead of being held at the prescribed value. Both now impose the full boundary velocity, matching Dirichlet #156
  • Fixed Slip contributing nothing to vector diffusion. No Laplacian{Linear} method was defined for vector fields, so the unqualified method returned zero coefficients for a slip patch. Slip now projects onto the face tangent plane and matches Symmetry coefficient for coefficient, and gains the missing Si method #156
  • Fixed Slip Divergence{BoundedUpwind} dropping the face-normal component for vector fields #156
  • Fixed simulation start-up scaling with the number of boundary patches instead of mesh size. Patch sizes were passed to the KernelAbstractions kernel constructor, where they become StaticSize type parameters, so every kernel launched once per patch was compiled again for each distinct patch size. On the motorBike case (354k cells, 72 patches) the first SIMPLE iteration spent 259 s in compilation, against 51 s for a 1.68M cell case with 6 patches, which presents as the solver hanging after Starting SIMPLE loops.... Sizes are now passed at launch through _dynamic_setup, so all patches share one compiled kernel and the same first iteration takes 73 s #158
  • Fixed wall functions racing on cells that own more than one wall face. set_production!, constrain! and the mixing-length correct_nut_wall! assigned the cell value directly, so with several faces writing the same cell the result depended on which face won, both within a patch and across patches, and varied from run to run. They now sum the contributions of every face that touches the cell and divide by their number, so each face carries equal weight and the result no longer depends on the order the faces are processed in. On the motorBike case 3220 of 45454 wall cells own 2 or 3 wall faces; single-face cells are unchanged and residuals are now reproducible run to run #158
  • Fixed the mass flux correction using the relaxed pressure field. The pressure equation is solved for the full predicted continuity error, so correcting the flux with a relaxed pressure removes only a fraction of it and leaves a residual imbalance that does not vanish as the solution converges. The flux is now corrected from the unrelaxed solution and relaxation is applied afterwards, to the momentum correction only #158
  • Fixed Wall and RotatingWall Divergence{BoundedUpwind} returning a zero source for vector fields. The bounded scheme subtracts ap from the diagonal, so a zero source implicitly convects the cell's own velocity through the wall face instead of the prescribed boundary velocity. Both now carry the boundary velocity on the right-hand side #158
  • Fixed set_production! evaluating the KWallFunction production against a hard-coded zero wall velocity, which overpredicts the shear at a moving or rotating wall. The wall velocity is now taken from Uf #158
  • Fixed the pressure flux at Extrapolated pressure patches not being propagated to the boundary mass flux. The patch's Laplacian contribution is assembled from the owner-cell pressure available at assembly time, so once the pressure equation is solved the boundary flux is stale. The correction is now deferred and applied for patches whose velocity boundary condition is pressure-adjustable #158
  • Fixed the non-orthogonal face correction being applied to the pressure equation but not to the mass flux, so continuity did not close on non-orthogonal meshes. The correction is now stored per face and subtracted from the mass flux #158

Breaking

  • assign now requires every mesh boundary to be assigned exactly once and throws an ArgumentError naming any missing or duplicated patch. Previously only the number of boundary conditions was checked, so an assignment that named one patch twice and omitted another was accepted, leaving a patch without a boundary condition #154

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.6.0 - 2026-07-17

Added

  • Added optional adaptive time stepping based on Courant number control (AdaptiveTimeStepping) #98
  • Added runtime calculation of Q-criterion #107
  • Added explicit Godunov-type compressible solver (godunov!) for Supersonic flows using Rusanov (local Lax-Friedrichs) and HLLC flux schemes with first and second order spatial reconstruction (MinMod, VanLeer, Superbee limiters) and adaptive CFL-based time stepping. This is the recommended solver for high-speed flows in the range 0.8 < M < 5 approximately #112
  • Added experimental incompressible steady MFR solver #114
  • Added runtime probe extraction #115
  • Added initialisation logic for new Multiphase solver user-level API #117
  • New Eulerian Thin Film model solver (2D only) #120
  • Added setField_Expression! to SetFields utility to initialise function-based fields. #124
  • Added VanLeer, upwind, and gradient interpolation schemes for scalar and vector face fields #124
  • Added multiphase solver (VOF model only) with two supporting functionality tests. #132
  • Added mixture model inside multiphase solver with supporting unit and functionality tests. #136
  • Extended initialise! API with function-based overloads for ScalarField and VectorField#135
  • Pressure-based compressible solvers have been extended to include SensibleEnthalpy and InternalEnergy formulations. The solver now correctly handles Compressible fluids, activating the transonic correction that allows shock capturing for moderate high-speed applications, typically in the range 0.8 < M < 2 #145

Fixed

  • Add implementation of Periodic boundaries to handle the implicit source term - fixes operation of models that use Si terms #95
  • Fixed implementation of implicit boundaries in #96 which where missing atomics #100
  • Fixed calculation of the residuals to use the relative residual norm, norm(b - Ax)/norm(b), the numerator in this expression was calculated incorrectly previously, giving a 1/sqrt(n) relation (where n is the number of cells in the grid). whilst the operation of the solvers remains the same, user may find that convergence criteria may need to be increased (specially for larger grids)#102. Move calculation to scale residuals in #134
  • UNV2: Fix calculation of cell volumes and centroid for boundary cells was incorrect and missing boundary face contributions (only for 2D UNV meshes)#106
  • Fixed implementation of k-omega LKE transition model and how wall distance field is calculated to ensure it is GPU compatible #109
  • Fixed mixed-precision mesh conversion to preserve user-selected integer and floating-point types #125
  • Fixed pressure boundary mass-flux correction #125
  • Fixed turbulent effective viscosity updates so turbulence models include eddy viscosity again, reverting a regression introduced in #120 #125
  • Add missing constructor for LES{Laminar} supporting qDNS simulations #144

Changed

  • Improved stability of Periodic boundaries by making the implementation fully implicit #96
  • 4x speed improvement for the method construct_periodic #97
  • +50x speed improvement for the method construct_periodic and also more robust algorithm used #100
  • Implementation to correct mass flux uses matrix coefficients directly for better stability when using periodic boundary conditions #100
  • New method to enforce matrix symmetry of scalar model equations when the only term is a laplacian #100
  • Change calculation of face interpolation weights to use face normal aligned weights, this is more physical than the current method using face-based distances (in preparation for formal support for non-orthogonality correction)#101
  • 20x improvement loading and processing 3D UNV mesh files #106
  • Updated discretisation, boundary condition, turbulence, solver, and preconditioner paths to avoid unintended Float64 promotion on Float32 meshes #125
  • FOAM reader has been updated for robustness and speed (~4x faster) #126
  • Rewrote pressure_force and viscous_force as backend-agnostic (CPU/GPU) kernels, changed their signature to (patch, model, config), and corrected the pressure force to scale by reference density only for incompressible (kinematic pressure) flows #129
  • Improved calculation of viscous force using wall eddy viscosity and perpendicular cell-to-face distance #130

Breaking

  • Pressure-based compressible solvers (for WeaklyCompressible and Compressible fluids) have been extended to work with SensibleEnthalpy and InternalEnergy formulations. Thus, to capture this change, the energy field has been renamed he instead of h #145

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.5.2 - 2025-01-16

Added

  • Initial support for mixed precision (UNV meshes only) #67
  • New solver for simulating conduction in solids #65
  • New LES turbulent kinetic energy one equation model (KEquation) #71
  • Surface tension model for fluids #72
  • High fidelity viscosity models for H2 and N2 #72
  • High fidelity thermal conductivity models for H2 and N2 #72
  • SetFields utility that allows to set a field to desired value within a box / circle / sphere #73
  • Helmholtz Energy equation of state and supporting framework for it for H2 and N2 #75
  • FieldAverage and FieldRMS of Vector and Scalar fields #78
  • Added RotatingWall velocity boundary condition #81
  • Added KOmegaSST turbulence model #82
  • Added capability to write out TensorFields to .vtk/vtu including Reynolds Stress Tensor #84
  • Added capablity to write out TensorFields to .vtk/vtu including Reynolds Stress Tensor #87
  • Extended DirichLetFunction to accept functions defining boundary condition for ScalarFields #89
  • Implemented CrankNicolson time scheme (second order implicit-explicit) #90
  • Implementation of a new AMG solver for symmetric matrices #141

Fixed

  • The UNV3D_mesh reader has been updated to ensure that the ordering of face nodes is determined in a more robust manner. This resolves some issues when loading a UNV mesh that is later used to store simulation results in the OpenFOAM format #64
  • In the construction of a Physics object, the boundary_map function returned a boundary_info struct which was incorrectly using an abstract type Integer. This resulted in a failure to convert a Physics object to the cpu and back to the gpu #67
  • Fixed calculation of interpolation weights for 2D UNV grids #68
  • Added method for boundary interpolation when using LUST and DirichletFunction#68
  • Fixed boundary_interpolation for DirichletFunction when a function is passed #79

Changed

  • The constructors for ScalarField and FaceScalarField now include a store_mesh keyword argument to request a reference of the mesh to be stored (default) or not (setting store_mesh=false). This can be used to not include references to the mesh for each field in VectorFields and TensorFields. This has improved compile times and decreased simulation times (particularly on the GPU - perhaps due to freeing registers used to carry unnecessary type information) #69
  • Internally, the calculation of interpolation weights and other geometric properties are calculated using the same function (defined in the Mesh module) #69
  • The default discretisation for laplacian terms uses the over-relaxed formulation by default. This will have no effect on orthogonal grids, but tends to be more robust in complex geometries at the expense of accuracy, which can be recovered by adding additional orthogonal correction loops (using the key word argument ncorrectors in the run! function) #73
  • Cleaned code for all solvers and improved stability of incompressible solver by removing the update of the mass flow based on the velocity field from the previous iteration. The mass flow is now corrected directly from the latest pressure solution #76
  • Improved calculation of wall distance for turbulence models to prevent negative square roots resulting in DomainErrors #139

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.5.1 - 2025-07-03

Added

  • New Empty boundary condition allowing 2D simulations with OpenFOAM 2D-compatible grids#63
  • Tests for the Smagorinsky LES model have been included #63

Fixed

  • No fixes included

Changed

  • Transient simulation results for VTK and VTU files use the format time_<iteration> #63

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.5.0 - 2025-06-28

Added

  • New boundary conditions Extrapolated and Zerogradient have been added. Both assign a zero gradient boundary condition, however, their implementation differs. Extrapolated assigns the zero gradient condition semi-implicitly (using the cell centre unknown and the cell centre value from the previous iteration). Zerogradient assigns the gradient in the boundary faces explicitly #61
  • The workgroup size for the CPU backend can now be automatically chosen when workgroup=AutoTune(). This uses a simple ceiling division internally where the number of elements in the kernel are divided by the number of available threads. This results in a small 10% performance gain #61

Fixed

  • Fixed the implementation for the calculation of the wall distance to work on GPUs #49
  • The new approach for handling user-provided boundary conditions now allows extracting the wall velocity specified by the user for the Wall boundary condition #61

Changed

  • In the calculation of wall function properties the user-provided wall velocity is now used, instead of hard-coded to no-slip #49
  • The functions set_production!, set_cell_value! and correct_nut_wall!, in RANS_functions.jl have been updated, removing conditional branch used in the generated calling them. Now these functions use multiple dispatch to allow specialising the wall function framework to ease the development of new wall functions #57
  • New NeumannFunction has been created, mirroring the DirichletFunction, providing a Neumann boundary condition defined with a user-provided struct (a generic framework accepting a function to set the gradient at the boundary is not yet available) #57
  • update_user_boundary function, extension has been reverted, overwriting the changes made to expose the ModelEquation type to it in #55
  • User-provided boundary conditions are no longer stored within fields, instead a NamedTuple is constructed using the function assign, which is passed to solvers using the Configuration struct.
  • Configuration setting provided by the user at the top-level API are now stored in predefined structs, instead of using NamedTuples, this change should put less pressure on the compiler #61
  • Internally, kernel launches have been updated to use a static NDrange. This resulted in a 20% speed improvement on GPU backends (only NVidia GPUs tested) #61

Breaking

  • The definition of Krylov solvers in the previous API used types exported directly from Krylov.jl. Now solvers are defined using instances of types defined in XCALibre.jl. As an example, previously the CG solver was defined using the type CgSolver now this solver is defined using the instance Cg() where the suffix "Solver" has been dropped. This applies to all previously available solver choices #60
  • The Green-Gauss method for calculating the gradient is now Gauss which is more descriptive than the previous name Orthogonal
  • The internals for handling user-provided boundary conditions have been updated in preparation for extending the code for handling multiple regions. Thus, the syntax for assigning boundary conditions has changed. The most noticeable change is the removal of the the @assign! macro, replaced by the function assign. See the documentation for details #61
  • The set_solver, set_hardware and set_runtime top-level functions have been replaced with SolverSetup, Hardware and Runtime structures. This allowed storing user-provided setup information in structs instead of NamedTuples to reduce the burden on the compiler #61
  • The implementation of the FixedTemperature boundary condition has been simplified, resulting in a breaking change. The definition of the boundary condition now requires users to provide an Enthalphy model. This interface should make it easier to extend the current implementation to new forms of the energy equation #61

Deprecated

  • No functions deprecated

Removed

  • The functions set_solver, set_runtime, set_hardware and the macro @assign! have been removed/replaced #61

Version v0.4.2 - 2025-04-02

Added

  • A very simple 2D block mesh generator has been added (not ready for general use as it needs to be documented)#41
  • Implementation of Wall boundary conditions specialised for ScalarField #45
  • Simulation results can now be written to VTK and OpenFOAM formats. The format can be selected using the output keyword argument in the run! function. The formats available are VTK() and OpenFOAM() #47

Fixed

  • Fixed the implementation for the calculation of the wall distance #45

Changed

  • In preparation for hybrid RANS/LES models, the object Turbulence is now passed to the TurbulenceModel object to allow for a more general call of turbulence!. This changes the implementation of turbulence! for all models slightly #46
  • Improved stability on meshes with warped faces by changing how face normals are calculated. XCALibre now uses an area-weighted face normal calculation based on the decomposition of faces into triangles #47
  • Removed VTK module and moved functionality to a new IOFormats module #47
  • Internally the function model2vtk has been replaced with save_output within all solvers. The specialisation for writing ModelPhysics models to file has also changed to save_output. The arguements pass to this function have also changed, the name of the file is not passed to the function, instead, the current runtime variable time is required #47
  • The Momentum object, part of the ModelPhysics object, now includes the FaceScalarField and FaveVectorField in addition to existing fields in preparation for non-uniform boundary definition when using the OpenFOAM output format #47

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.4.1 - 2025-03-06

Added

  • No new functionality has been added

Fixed

  • Remove fragile precompile statements of the form var"#*#*" causing errors in Julia v1.12 #40

Changed

  • In preparation for implementation of hybrid models e.g. DES, the signature of the turbulence! function has been updated to use the AbstractTurbulenceModel #39
  • New semi-implicit implementation for periodic boundary conditions (implicit treatment of laplacian terms)#40

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.4.0 - 2025-02-17

Added

  • Implementation of Symmetry boundary condition for ScalarField types #30
  • New macro to help define boundary conditions that will dispatch to Scalar or VectorField types #30
  • Added eltype method for both Scalar and VectorField types to simplify the development of new kernels where type information is needed #30
  • New gradient limiters FaceBased and MFaceBased for limiting gradients based on cell faces, where MFaceBased is a multidimensional version, and it is generally recommended over FaceBased #30
  • Support for INTEL hardware #32

Fixed

  • Calling JacobiSmoother now works on the GPU #30
  • Implemented SparseXCSR as wrapper for SparseMatrixCSR on the CPU to resolve display/print errors #30
  • The convergence criteria for solvers is used consistently to stop/control simulation runtime #36
  • Consistent display of residuals on screen during simulations #36
  • Fixed calculation of wall distance and apply_boundary_conditions! arguments for LKE model #36

Changed

  • The calculation of gradients has been improved by merging computations into a single kernel, improving performance of gradient kernels by around 10-30%, most noticable for vector gradients #30
  • Improved calculation of non-orthogonal calculation (more tests are still needed), although tests have proven to be stable #30
  • Improved documentation/readme on supported GPU backends/hardware and make users aware of potential F32 limitation on some hardware

Breaking

  • The top level API for all solvers no longer takes the keyword argument limit_gradient for activating gradient limiter. New gradient limiters have been added and can be selected/configured when assigning numerical schemes with the Schemes function #30

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.3.3 - 2024-12-24

Added

  • Added experimental support for NVIDIA ILU0 and IC0 preconditioners #23
  • Added JacobiSmoother that can be used with all linear solvers (improving initial guess) #23
  • New function activate_multithread is available to set up matrix-vector multiplication in parallel on the CPU
  • Initial benchmark added in the documentation

Fixed

  • No fixes

Changed

  • Internally the sparse matrix format has been changed to CSR. This has improved performance by 2x-4x (case dependent)
  • Multithreaded sparse matrix vector multiplication is now functional #23
  • Precompilation errors on Julia v1.11 addressed by bringing code from ThreadedSparseCSR.jl #24
  • Update compat entry for Atomix.jl to v1.0,0 #24
  • DILU preconditioner is now implemented to work with sparse matrices in CSR format and uses a hybrid approach (running on CPU) to allow working when using GPU backends (further work needed) #26
  • The implementation of RANS models and wall functions has been improved for consistency (resulting in some computational gains). The calculation of yPlusLam is only done once when constructing the wall function objects. The calculation of the velocity gradient is now only done within the turbulence model main function (turbulence!) #28

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • The ILU0 and LDL preconditioners has been temporarily removed #23

Version v0.3.2 - 2024-11-08

Added

  • Added support for CI testing, Dependabot, and CompatHelper

Fixed

  • Fixed tests for mesh conversion and standardised tolerances for test checks of incompressible solvers #16

Changed

  • No changes

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.3.1 - 2024-10-18

Added

  • Vastly improved documentation with new examples provided #12
  • Changelog added to record changes more clearly. Record kept in Release notes

Fixed

  • The calculation of gradients can be limited for stability. This functionality can be activated by passing the key work argument limit_gradient to the run! function. The implementation has been improved for robustness #12
  • Removed face information being printed when Mesh objects are created to stop printing a ERROR: Scalar indexing is disallowed message #13

Changed

  • Master branch protected and requires PRs to push changes

Breaking

  • No breaking changes

Deprecated

  • No functions deprecated

Removed

  • No functionality has been removed

Version v0.3.0 - 2024-09-21

  • New name - XCALibre.jl - which is now registered in the General Julia registry
  • Can do 3D and GPU accelerated simulations
  • Can read .unv and OpenFOAM mesh files (3D)
  • Can do incompressible and compressible simulations
  • RANS and LES models available
  • User-provided functions or neural networks for boundary conditions
  • Reasonably complete "user" documentation now provided
  • Made repository public (in v0.2 the work was kept in a private repository and could only do 2D simulations)
  • Tidy up mesh type definitions by @mberto79 in #5
  • Adapt code base to work with new mesh format by @mberto79 in #6
  • Mesh boundary struct changes PR by @TomMazin in #7
  • Mesh boundary struct changes PR fix by @TomMazin in #8

Version v0.2.0 - 2023-01-23

  • New mesh format and type implemented that are GPU friendly.
  • No functionality changes

Version v0.1.0 - 2023-01-23

Initial release

2D implementation of classic incompressible solvers for laminar and turbulent flows:

  • Framework for equation definition
  • SIMPLE nad PISO algorithms
  • Read UNV meshes in 2D
  • Capability for RANS models
  • Various discretisation schemes available
  • Planned extension to 3D and GPU acceleration!