Merge pull request #102 from mberto79/HM/Fix-L2-norm-residual-calculation

Fix residual calculation

Author: mberto79

CHANGELOG.md

@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### 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)
+* 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](@ref)
 
 ### Changed
 * Improved stability of `Periodic` boundaries by making the implementation fully implicit [#96](@ref)

examples/2D_cylinder_heated.jl

@@ -20,7 +20,7 @@ mesh_dev = adapt(backend, mesh)
 
 velocity = [0.5, 0.0, 0.0]
 noSlip = [0.0, 0.0, 0.0]
-nu = 1e-2
+nu = 1e-4
 Re = (0.2*velocity[1])/nu
 gamma = 1.4
 cp = 1005.0
@@ -46,24 +46,24 @@ BCs = assign(
     (
         U = [
             Dirichlet(:inlet, velocity),
-            Neumann(:outlet, 0.0),
+            Zerogradient(:outlet),
             Wall(:cylinder, noSlip),
             Symmetry(:bottom, 0.0),
             Symmetry(:top, 0.0)
         ],
         p = [
-            Neumann(:inlet, 0.0),
+            Zerogradient(:inlet),
             Dirichlet(:outlet, pressure),
             Wall(:cylinder, 0.0),
-            Neumann(:bottom, 0.0),
-            Neumann(:top, 0.0)
+            Zerogradient(:bottom),
+            Zerogradient(:top)
         ],
         h = [
             FixedTemperature(:inlet, T=300.0, Enthalpy(cp=cp, Tref=288.15)),
-            Neumann(:outlet, 0.0),
+            Zerogradient(:outlet),
             FixedTemperature(:cylinder, T=330.0, Enthalpy(cp=cp, Tref=288.15)),
-            Neumann(:bottom, 0.0),
-            Neumann(:top, 0.0)
+            Zerogradient(:bottom),
+            Zerogradient(:top)
         ]
     )
 )

examples/2D_cylinder_heated_unsteady.jl

@@ -45,24 +45,24 @@ BCs =assign(
     (
         U = [
             Dirichlet(:inlet, velocity),
-            Neumann(:outlet, 0.0),
+            Zerogradient(:outlet, 0.0),
             Wall(:cylinder, noSlip),
-            Neumann(:bottom, 0.0),
-            Neumann(:top, 0.0)
+            Zerogradient(:bottom, 0.0),
+            Zerogradient(:top, 0.0)
         ],
         p = [
-            Neumann(:inlet, 0.0),
+            Zerogradient(:inlet, 0.0),
             Dirichlet(:outlet, pressure),
-            Neumann(:cylinder, 0.0),
-            Neumann(:bottom, 0.0),
-            Neumann(:top, 0.0)
+            Zerogradient(:cylinder, 0.0),
+            Zerogradient(:bottom, 0.0),
+            Zerogradient(:top, 0.0)
         ],
         h = [
             FixedTemperature(:inlet, T=temp, Enthalpy(cp=cp, Tref=288.15)),
-            Neumann(:outlet, 0.0),
+            Zerogradient(:outlet, 0.0),
             FixedTemperature(:cylinder, T=330.0, Enthalpy(cp=cp, Tref=288.15)),
-            Neumann(:bottom, 0.0),
-            Neumann(:top, 0.0)
+            Zerogradient(:bottom, 0.0),
+            Zerogradient(:top, 0.0)
         ]
     )
 )

examples/2D_flatplate_lowRe.jl

@@ -35,8 +35,8 @@ BCs = assign(
             Dirichlet(:inlet, velocity),
             Zerogradient(:outlet),
             Wall(:wall, [0.0, 0.0, 0.0]),
-            Extrapolated(:top)
-            # Zerogradient(:top)
+            # Extrapolated(:top)
+            Zerogradient(:top)
             # Symmetry(:top)
         ],
         p = [
@@ -51,16 +51,16 @@ BCs = assign(
             Dirichlet(:inlet, k_inlet),
             Zerogradient(:outlet),
             Dirichlet(:wall, 0.0),
-            Extrapolated(:top)
-            # Zerogradient(:top)
+            # Extrapolated(:top)
+            Zerogradient(:top)
             # Symmetry(:top)
         ],
         omega = [
             Dirichlet(:inlet, ω_inlet),
             Zerogradient(:outlet),
             OmegaWallFunction(:wall),
-            Extrapolated(:top)
-            # Zerogradient(:top)
+            # Extrapolated(:top)
+            Zerogradient(:top)
             # Symmetry(:top)
         ],
         nut = [
@@ -78,7 +78,8 @@ schemes = (
     U = Schemes(divergence=Upwind),
     p = Schemes(divergence=Upwind),
     k = Schemes(divergence=Upwind),
-    omega = Schemes(divergence=Upwind)
+    omega = Schemes(divergence=Upwind),
+    y = Schemes()
 )
 
 
@@ -87,29 +88,36 @@ solvers = (
         solver      = Bicgstab(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(), 
         convergence = 1e-7,
-        relax       = 0.7,
+        relax       = 0.8,
     ),
     p = SolverSetup(
         solver      = Cg(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(), 
         convergence = 1e-7,
         relax       = 0.2,
     ),
+    y = SolverSetup(
+        solver      = Cg(), # Bicgstab(), Gmres()
+        preconditioner = Jacobi(), 
+        convergence = 1e-7,
+        relax       = 0.9,
+        itmax = 5000
+    ),
     k = SolverSetup(
         solver      = Bicgstab(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(), 
         convergence = 1e-7,
-        relax       = 0.7,
+        relax       = 0.6,
     ),
     omega = SolverSetup(
         solver      = Bicgstab(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(), 
         convergence = 1e-7,
-        relax       = 0.7,
+        relax       = 0.6,
     )
 )
 
-runtime = Runtime(iterations=1000, write_interval=100, time_step=1)
+runtime = Runtime(iterations=5000, write_interval=100, time_step=1)
 
 config = Configuration(
     solvers=solvers, schemes=schemes, runtime=runtime, hardware=hardware, boundaries=BCs)

examples/3D_cascade_periodic.jl

@@ -76,25 +76,25 @@ solvers = (
     U = SolverSetup(
         solver      = Bicgstab(), #Cg(), # Bicgstab(), Gmres(), #Cg()
         preconditioner = Jacobi(),
-        convergence = 1e-8,
+        convergence = 1e-7,
         # # transient setup
         # atol = 1e-6,
         # relax=1
 
         # steady setup
-        relax       = 0.6,
+        relax       = 0.7,
         rtol = 1e-3
     ),
     p = SolverSetup(
         solver      = Cg(), #Gmres(), #Cg(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(),
-        convergence = 1e-8,
+        convergence = 1e-7,
         # # transient setup
         # atol = 1e-6,
         # relax=1
 
         # steady setup
-        relax       = 0.15,
+        relax       = 0.3,
         rtol = 1e-3
     )
 )
@@ -113,9 +113,9 @@ initialise!(model.momentum.p, 0.0)
 
 residuals = run!(model, config, output=OpenFOAM()) # 353 iterations!
 
-# using Plots
-# fig = plot(; xlims=(0,runtime.iterations), ylims=(1e-10, 1e-4))
-# plot!(fig, 1:runtime.iterations, residuals.Ux, yscale=:log10, label="Ux")
-# plot!(fig, 1:runtime.iterations, residuals.Uy, yscale=:log10, label="Uy")
-# plot!(fig, 1:runtime.iterations, residuals.p, yscale=:log10, label="p")
-# fig
+using Plots
+fig = plot(; xlims=(0,runtime.iterations), ylims=(1e-10, 1e-1))
+plot!(fig, 1:runtime.iterations, residuals.Ux, yscale=:log10, label="Ux")
+plot!(fig, 1:runtime.iterations, residuals.Uy, yscale=:log10, label="Uy")
+plot!(fig, 1:runtime.iterations, residuals.p, yscale=:log10, label="p")
+fig

src/Solve/Solve_1_api.jl

@@ -406,13 +406,12 @@ function residual(eqn, component, config)
 
     # Previous definition
     Fx .= A * values
-    # @inbounds @. R = (b - Fx)^2
     xcal_foreach(R, config) do i 
-            R[i] = (b[i] - Fx[i])^2
+            @inbounds R[i] = (b[i] - Fx[i])^2
     end
     normb = norm(b)
-    denominator = ifelse(normb>0,normb, 1)
-    Residual = sqrt(mean(R)) / denominator
+    denominator = ifelse(normb > eps(normb), normb, one(normb))
+    Residual = sqrt(sum(R)) / denominator
     return Residual
 end