Merge branch 'main' into HM/normal-aligned-interpolation-weights

Author: mberto79

CHANGELOG.md

@@ -3,18 +3,22 @@
 The format used for this `changelog` is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). 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. 
 
-## Version [v0.5.2] - 2025-01-16
+## Version [v0.5.3] - 2026-02-19
 
 ### Added
 *  No new functionality added
 
 
 ### Fixed
 * Add implementation of `Periodic` boundaries to handle the implicit source term - fixes operation of models that use `Si` terms [#95](@ref)
+* Fixed implementation of implicit boundaries in [#96](@ref) which where missing atomics [#100](@ref)
 
 ### Changed
 * Improved stability of `Periodic` boundaries by making the implementation fully implicit [#96](@ref)
 * 4x speed improvement for the method `construct_periodic` [#97](@ref)
+* +50x speed improvement for the method `construct_periodic` and also more robust algorithm used [#100](@ref)
+* Implementation to correct mass flux uses matrix coefficients directly for better stability when using periodic boundary conditions [#100](@ref)
+* New method to enforce matrix symmetry of scalar model equations when the only term is a laplacian [#100](@ref)
 
 ### Breaking
 * No breaking changes

src/Solve/Solve_1_api.jl

@@ -169,7 +169,9 @@ function solve_equation!(
 
     discretise!(eqn, phi, config)       
     apply_boundary_conditions!(eqn, phiBCs, nothing, time, config)
-    make_symmetric!(eqn, config) # added this to test stability of periodic boundaries
+    if length(eqn.model.terms) == 1 && typeof(eqn.model.terms[1]) <: Laplacian
+        make_symmetric!(eqn, config) # added this to test stability of periodic boundaries
+    end
     setReference!(eqn, ref, 1, config)
     if !isnothing(irelax)
         implicit_relaxation!(eqn, phi.values, irelax, nothing, config)

src/Solvers/Solvers_1_CSIMPLE.jl

@@ -226,6 +226,7 @@ function CSIMPLE(
         # Pressure correction
         inverse_diagonal!(rD, U_eqn, config)
         interpolate!(rhorDf, rD, config)
+        correct_interpolation_periodic(rhorDf, rD, boundaries.U, config)
         @. rhorDf.values *= rhof.values
 
         remove_pressure_source!(U_eqn, ∇p, config)
@@ -297,11 +298,11 @@ function CSIMPLE(
 
         if typeof(model.fluid) <: Compressible
             # @. mdotf.values += (pconv.values*(pf.values) - pgrad.values*rhorDf.values)  
-            correct_mass_flux(mdotf, p, rhorDf, config)
+            correct_mass_flux(mdotf, p_eqn, config)
             @. mdotf.values += pconv.values*(pf.values)
         elseif typeof(model.fluid) <: WeaklyCompressible
             # @. mdotf.values -= pgrad.values*rhorDf.values
-            correct_mass_flux(mdotf, p, rhorDf, config)
+            correct_mass_flux(mdotf, p_eqn, config)
         end
 
         correct_velocity!(U, Hv, ∇p, rD, config)

src/Solvers/Solvers_1_SIMPLE.jl

@@ -75,17 +75,10 @@ function setup_incompressible_solvers(
 
     p_eqn = (
         - Laplacian{schemes.p.laplacian}(rDf, p) == - Source(divHv)
-        # Laplacian{schemes.p.laplacian}(rDf, p) == Source(divHv)
-        # Laplacian{schemes.p.laplacian}(rDf, p) == Source(divHv)
     ) → ScalarEquation(p, boundaries.p)
 
     @info "Initialising preconditioners..."
 
-    # @reset U_eqn.preconditioner = set_preconditioner(
-    #                 solvers.U.preconditioner, U_eqn, boundaries.U, config)
-    # @reset p_eqn.preconditioner = set_preconditioner(
-    #                 solvers.p.preconditioner, p_eqn, boundaries.p, config)
-
     @reset U_eqn.preconditioner = set_preconditioner(solvers.U.preconditioner, U_eqn)
     @reset p_eqn.preconditioner = set_preconditioner(solvers.p.preconditioner, p_eqn)
 
@@ -141,7 +134,6 @@ function SIMPLE(
 
     # Pre-allocate auxiliary variables
     TF = _get_float(mesh)
-    # prev = _convert_array!(prev, backend) 
     prev = KernelAbstractions.zeros(backend, TF, n_cells) 
 
     # Pre-allocate vectors to hold residuals 
@@ -212,17 +204,7 @@ function SIMPLE(
             limit_gradient!(schemes.p.limiter, ∇p, p, config)
         end
 
-        # explicit_relaxation!(p, prev, solvers.p.relax, config)
-
-        # Velocity and boundaries correction
-
-        # old approach
-        # correct_velocity!(U, Hv, ∇p, rD, config)
-        # interpolate!(Uf, U, config)
-        # correct_boundaries!(Uf, U, boundaries.U, time, config)
-        # flux!(mdotf, Uf, config) 
-
-        # new approach
+        # correct mass flux and velocity
         correct_mass_flux(mdotf, p_eqn, config)
         correct_velocity!(U, Hv, ∇p, rD, config)
 
@@ -278,7 +260,7 @@ function SIMPLE(
     return (Ux=R_ux, Uy=R_uy, Uz=R_uz, p=R_p)
 end
 
-### TEMP LOCATION FOR PROTOTYPING - NONORTHOGONAL CORRECTION 
+### TEMP LOCATION FOR PROTOTYPING
 
 function nonorthogonal_face_correction(eqn, grad, flux, config)
     mesh = grad.mesh
@@ -331,19 +313,11 @@ end
     T_hat = Sf - Ef # original
     faceCorrection = flux[fID]*gradf⋅T_hat
 
-    Atomix.@atomic b[cID1] += faceCorrection #*cell1.volume
-    Atomix.@atomic b[cID2] -= faceCorrection #*cell2.volume # should this be -ve?
-
-    # Atomix.@atomic b[cID1] -= faceCorrection #*cell1.volume
-    # Atomix.@atomic b[cID2] += faceCorrection #*cell2.volume # should this be -ve?
-        
+    Atomix.@atomic b[cID1] += faceCorrection
+    Atomix.@atomic b[cID2] -= faceCorrection 
+      
 end
 
-# +- => good match
-# -+ => looks worse at edges for gradient
-# -- => looks bad on top-right corner for gradient
-# ++ => looks bad on left grad and oscillations on the right
-
 function correction_weight(cells, faces, fi)
     (; ownerCells, centre) = faces[fi]
     cID1 = ownerCells[1]
@@ -405,31 +379,12 @@ end
         # need to get aN from sparse system
         zID = spindex(rowptr, colval, cID1, cID2)
         aN = nzval[zID]
-        # mdotf[fID] += -aN*(p2 - p1) # this worked well with periodics
         mdotf[fID] += aN*(p2 - p1) # positive because pressure eqn has negative sign
     end
 end
 
 ### Correct mass flux at periodic boundaries
 
-# function correct_mass_periodic(mdotf, p_eqn, pBCs, config)
-#     (; faces, cells, boundary_cellsID) = mdotf.mesh
-#     (; hardware) = config
-#     (; backend, workgroup) = hardware
-
-#     p = p_eqn.model.terms[1].phi
-#     A = _A(p_eqn)
-#     nzval = _nzval(A)
-#     colval = _colval(A)
-#     rowptr = _rowptr(A)
-
-#     for BC ∈ pBCs
-#         _correct_mass_periodic_dispatch(
-#             BC, mdotf, p, nzval, colval, rowptr, cells, faces, backend, workgroup)
-#     end
-
-# end
-
 correct_mass_periodic(arg...) = nothing
 
 function correct_mass_periodic(
@@ -463,7 +418,6 @@ end
 
 end
 
-
 ### Correct interpolation at periodic boundaries
 
 function correct_interpolation_periodic(phif, phi, BCs, config)

src/Solvers/Solvers_2_CPISO.jl

@@ -249,6 +249,7 @@ function CPISO(
         # Pressure correction
         inverse_diagonal!(rD, U_eqn, config)
         interpolate!(rhorDf, rD, config)
+        correct_interpolation_periodic(rhorDf, rD, boundaries.U, config)
         @. rhorDf.values *= rhof.values
 
         remove_pressure_source!(U_eqn, ∇p, config)
@@ -313,10 +314,10 @@ function CPISO(
             end
 
             if typeof(model.fluid) <: Compressible
-                correct_mass_flux(mdotf, p, rhorDf, config)
+                correct_mass_flux(mdotf, p_eqn, config)
                 @. mdotf.values += pconv.values*(pf.values)
             elseif typeof(model.fluid) <: WeaklyCompressible
-                correct_mass_flux(mdotf, p, rhorDf, config)
+                correct_mass_flux(mdotf, p_eqn, config)
             end
 
             # TO-DO: this needs to be exposed to users eventually
@@ -325,10 +326,6 @@ function CPISO(
 
             # Velocity and boundaries correction
             correct_velocity!(U, Hv, ∇p, rD, config) # why is this not rhorD?
-
-            # Lines below should not be needed, interpolation to Uf happens in grad calcs
-            # interpolate!(Uf, U, config)
-            # correct_boundaries!(Uf, U, boundaries.U, time, config)
 
             @. dpdt.values = (p.values-prev)/runtime.dt
 

src/Solvers/Solvers_2_PISO.jl

@@ -161,12 +161,6 @@ function PISO(
                 limit_gradient!(schemes.p.limiter, ∇p, p, config)
             end
 
-            # old approach - keep for now!
-            # correct_velocity!(U, Hv, ∇p, rD, config)
-            # interpolate!(Uf, U, config)
-            # correct_boundaries!(Uf, U, boundaries.U, time, config)
-            # flux!(mdotf, Uf, config) # old approach
-
             # new approach
             correct_mass_flux(mdotf, p_eqn, config)
             correct_velocity!(U, Hv, ∇p, rD, config)