compensated summation for RK4, Diff and tracer

Author: milankl

src/constants.jl

@@ -18,10 +18,10 @@ function SSPRK3coeff{T}(P::Parameter,Δt_Δ::T) where T
     s = n^2
     kn = n*(n+1) ÷ 2 + 1
     mn = (n-1)*(n-2) ÷ 2 + 1
-    Δt_Δn = T(Δt_Δ/(n^2-n))
-    kna = T((n-1)/(2n-1))
-    knb = T(n/(2n-1))
-    Δt_Δnc = T(Δt_Δ/(n*(2n-1)))
+    Δt_Δn = convert(T,Δt_Δ/(n^2-n))
+    kna = convert(T,(n-1)/(2n-1))
+    knb = convert(T,n/(2n-1))
+    Δt_Δnc = convert(T,Δt_Δ/(n*(2n-1)))
 
     return SSPRK3coeff{T}(n,s,kn,mn,Δt_Δn,kna,knb,Δt_Δnc)
 end
@@ -76,37 +76,41 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
     end
 
     # Δt/(s-1) for SSPRK2
-    Δt_Δs = Tprog(G.dtint/G.Δ/(P.RKs-1))
+    Δt_Δs = convert(Tprog,G.dtint/G.Δ/(P.RKs-1))
 
     # time step and half the time step including the grid spacing as this is not included in the RHS
-    Δt_Δ = Tprog(G.dtint/G.Δ)
-    Δt_Δ_half = Tprog(G.dtint/G.Δ/2)
+    Δt_Δ = convert(Tprog,G.dtint/G.Δ)
+    Δt_Δ_half = convert(Tprog,G.dtint/G.Δ/2)
 
     # coefficients for SSPRK3
     SSPRK3c = SSPRK3coeff{Tprog}(P,Δt_Δ)
 
     # BOUNDARY CONDITIONS AND PHYSICS
-    one_minus_α = Tprog(1-P.α)      # for the ghost point copy/tangential boundary conditions
-    g = T(P.g)                      # gravity - for Bernoulli potential
+    one_minus_α = convert(Tprog,1-P.α)      # for the ghost point copy/tangential boundary conditions
+    g = convert(T,P.g)                      # gravity - for Bernoulli potential
 
-    # BOTTOM FRICTION COEFFICENTS
+    # BOTTOM FRICTION COEFFICIENTS
     # incl grid spacing Δ for non-dimensional gradients
     # include scale for quadratic cD only to unscale the scale^2 in u^2
-    cD = T(-G.Δ*P.cD/P.scale)     # quadratic drag [m]
-    rD = T(-G.Δ/(P.τD*24*3600))   # linear drag [m/s]
+    cD = convert(T,-G.Δ*P.cD/P.scale)     # quadratic drag [m]
+    rD = convert(T,-G.Δ/(P.τD*24*3600))   # linear drag [m/s]
 
     # INTERFACE RELAXATION FREQUENCY
     # incl grid spacing Δ for non-dimensional gradients
-    γ = T(G.Δ/(P.t_relax*3600*24))    # [m/s]
+    γ = convert(T,G.Δ/(P.t_relax*3600*24))    # [m/s]
 
     # BIHARMONIC DIFFUSION
     # undo scaling here as smagorinksy diffusion contains scale^2 due to ~u^2
-    cSmag = T(-P.cSmag/P.scale)   # Smagorinsky coefficient
-    νB = T(-P.νB/30000)           # linear scaling based on 540m^s/s at Δ=30km
+    cSmag = convert(T,-P.cSmag/P.scale)   # Smagorinsky coefficient
+    νB = convert(T,-P.νB/30000)           # linear scaling based on 540m^s/s at Δ=30km
 
     # TRACER ADVECTION
-    τSST = T(G.dtadvint/(P.τSST*3600*24))    # tracer restoring [1]
-    jSST = T(G.dtadvint/(P.jSST*3600*24))    # tracer consumption [1]
+    τSST = convert(T,G.dtadvint/(P.τSST*3600*24))   # tracer restoring [1]
+    jSST = convert(T,G.dtadvint/(P.jSST*3600*24))   # tracer consumption [1]
+
+    @unpack tracer_relaxation, tracer_consumption = P
+    τSST = tracer_relaxation ? τSST : zero(T)       # set zero as τ,j will be added   
+    jSST = tracer_consumption ? jSST : zero(T)      # and executed in one loop
 
     # TIME DEPENDENT FORCING
     ωFη = -2π*P.ωFη/24/365.25/3600
@@ -123,4 +127,4 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
                                 g,cD,rD,γ,cSmag,νB,τSST,jSST,
                                 ωFη,ωFx,ωFy,
                                 scale,scale_inv,scale_sst)
-end
+end

src/default_parameters.jl

@@ -53,16 +53,17 @@
     wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing
 
     # TIME STEPPING OPTIONS
-    time_scheme::String="SSPRK3"        # Runge-Kutta ("RK") or strong-stability preserving RK
+    time_scheme::String="RK"            # Runge-Kutta ("RK") or strong-stability preserving RK
                                         # "SSPRK2","SSPRK3","4SSPRK3"
     RKo::Int=4                          # Order of the RK time stepping scheme (2, 3 or 4)
     RKs::Int=3                          # Number of stages for SSPRK2
     RKn::Int=5                          # n^2 = s = Number of stages  for SSPRK3
-    cfl::Real=4.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
+    cfl::Real=0.9                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
     Ndays::Real=200.0                   # number of days to integrate for
     nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
     nstep_advcor::Int=0                 # advection and coriolis update every nstep_advcor time steps.
                                         # 0 means it is included in every RK4 substep
+    compensated::Bool=false             # Compensated summation in the time integration?
 
     # BOUNDARY CONDITION OPTIONS
     bc::String="periodic"               # "periodic" or anything else for nonperiodic

src/diffusion.jl

@@ -218,23 +218,36 @@ function viscous_tensor_constant!(  Diag::DiagnosticVars,
 end
 
 """Update u with bottom friction tendency (Bu,Bv) and biharmonic viscosity."""
-function add_drag_diff_tendencies!( u::Array{Tprog,2},
-                                    v::Array{Tprog,2},
+function add_drag_diff_tendencies!( u::Matrix{Tprog},
+                                    v::Matrix{Tprog},
                                     Diag::DiagnosticVars{T,Tprog},
                                     S::ModelSetup{T,Tprog}) where {T,Tprog}
 
     @unpack Bu,Bv = Diag.Bottomdrag
     @unpack LLu1,LLu2,LLv1,LLv2 = Diag.Smagorinsky
     @unpack halo,ep,Δt_diff = S.grid
+    @unpack compensated = S.parameters
+    @unpack du_comp,dv_comp = Diag.Tendencies
 
     m,n = size(u) .- (2*halo,2*halo)
     @boundscheck (m+2-ep,n+2) == size(Bu) || throw(BoundsError())
     @boundscheck (m,n+2) == size(LLu1) || throw(BoundsError())
     @boundscheck (m+2-ep,n) == size(LLu2) || throw(BoundsError())
 
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            u[i+2,j+2] += Δt_diff*(Tprog(Bu[i+1-ep,j+1]) + Tprog(LLu1[i,j+1]) + Tprog(LLu2[i+1-ep,j]))
+    if compensated
+        @inbounds for j ∈ 1:n
+            for i ∈ 1:m
+                du = Δt_diff*convert(Tprog,Bu[i+1-ep,j+1] + LLu1[i,j+1] + LLu2[i+1-ep,j]) - du_comp[i+2,j+2]
+                u_new = u[i+2,j+2] + du
+                du_comp[i+2,j+2] = (u_new - u[i+2,j+2]) - du
+                u[i+2,j+2] = u_new
+            end
+        end
+    else
+        @inbounds for j ∈ 1:n
+            for i ∈ 1:m
+                u[i+2,j+2] += Δt_diff*(Tprog(Bu[i+1-ep,j+1]) + Tprog(LLu1[i,j+1]) + Tprog(LLu2[i+1-ep,j]))
+            end
         end
     end
 
@@ -243,9 +256,20 @@ function add_drag_diff_tendencies!( u::Array{Tprog,2},
     @boundscheck (m,n+2) == size(LLv1) || throw(BoundsError())
     @boundscheck (m+2,n) == size(LLv2) || throw(BoundsError())
 
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-             v[i+2,j+2] += Δt_diff*(Tprog(Bv[i+1,j+1]) + Tprog(LLv1[i,j+1]) + Tprog(LLv2[i+1,j]))
+    if compensated
+        @inbounds for j ∈ 1:n
+            for i ∈ 1:m 
+                dv = Δt_diff*convert(Tprog,Bv[i+1,j+1] + LLv1[i,j+1] + LLv2[i+1,j]) - dv_comp[i+2,j+2]
+                v_new = v[i+2,j+2] + dv
+                dv_comp[i+2,j+2] = (v_new - v[i+2,j+2]) - dv
+                v[i+2,j+2] = v_new
+            end
+        end
+    else
+        @inbounds for j ∈ 1:n
+            for i ∈ 1:m
+                v[i+2,j+2] += Δt_diff*(Tprog(Bv[i+1,j+1]) + Tprog(LLv1[i,j+1]) + Tprog(LLv2[i+1,j]))
+            end
         end
     end
 end

src/preallocate.jl

@@ -60,6 +60,16 @@ end
     du::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)     # tendency of u without time step
     dv::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)     # tendency of v without time step
     dη::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)     # tendency of η without time step
+
+    # sum of tendencies (incl time step) over all sub-steps
+    du_sum::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo) 
+    dv_sum::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)
+    dη_sum::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)
+
+    # compensation for tendencies (variant of Kahan summation)
+    du_comp::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo) 
+    dv_comp::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)
+    dη_comp::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)
 end
 
 """Generator function for Tendencies VarCollection."""
@@ -422,6 +432,9 @@ end
 
     ssti::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst interpolated on departure points
     sst_ref::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst initial conditions for relaxation
+
+    # compensated summation
+    dsst_comp::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty)
 end
 
 """Generator function for SemiLagrange VarCollection."""

src/time_integration.jl

@@ -7,10 +7,13 @@ function time_integration(  Prog::PrognosticVars{Tprog},
     @unpack u0,v0,η0 = Diag.RungeKutta
     @unpack u1,v1,η1 = Diag.RungeKutta
     @unpack du,dv,dη = Diag.Tendencies
+    @unpack du_sum,dv_sum,dη_sum = Diag.Tendencies
+    @unpack du_comp,dv_comp,dη_comp = Diag.Tendencies
+
     @unpack um,vm = Diag.SemiLagrange
 
     @unpack dynamics,RKo,RKs,tracer_advection = S.parameters
-    @unpack time_scheme = S.parameters
+    @unpack time_scheme,compensated = S.parameters
     @unpack RKaΔt,RKbΔt = S.constants
     @unpack Δt_Δ,Δt_Δs = S.constants
 
@@ -55,6 +58,12 @@ function time_integration(  Prog::PrognosticVars{Tprog},
 
         if time_scheme == "RK"   # classic RK4,3 or 2
 
+            if compensated
+                fill!(du_sum,zero(Tprog))
+                fill!(dv_sum,zero(Tprog))
+                fill!(dη_sum,zero(Tprog))
+            end
+
             for rki = 1:RKo
                 if rki > 1
                     ghost_points!(u1,v1,η1,S)
@@ -74,10 +83,30 @@ function time_integration(  Prog::PrognosticVars{Tprog},
                     caxb!(η1,η,RKbΔt[rki],dη)   #η1 .= η .+ RKb[rki]*Δt*dη
                 end
 
-                # sum RK-substeps on the go
-                axb!(u0,RKaΔt[rki],du)          #u0 .+= RKa[rki]*Δt*du
-                axb!(v0,RKaΔt[rki],dv)          #v0 .+= RKa[rki]*Δt*dv
-                axb!(η0,RKaΔt[rki],dη)          #η0 .+= RKa[rki]*Δt*dη
+                if compensated      # accumulate tendencies
+                    axb!(du_sum,RKaΔt[rki],du)   
+                    axb!(dv_sum,RKaΔt[rki],dv)
+                    axb!(dη_sum,RKaΔt[rki],dη)
+                else    # sum RK-substeps on the go
+                    axb!(u0,RKaΔt[rki],du)          #u0 .+= RKa[rki]*Δt*du
+                    axb!(v0,RKaΔt[rki],dv)          #v0 .+= RKa[rki]*Δt*dv
+                    axb!(η0,RKaΔt[rki],dη)          #η0 .+= RKa[rki]*Δt*dη
+                end
+            end
+
+            if compensated
+                # add compensation term to total tendency
+                axb!(du_sum,-1,du_comp)             
+                axb!(dv_sum,-1,dv_comp)
+                axb!(dη_sum,-1,dη_comp)
+
+                axb!(u0,1,du_sum)   # update prognostic variable with total tendency
+                axb!(v0,1,dv_sum)
+                axb!(η0,1,dη_sum)
+                
+                dambmc!(du_comp,u0,u,du_sum)    # compute new compensation
+                dambmc!(dv_comp,v0,v,dv_sum)
+                dambmc!(dη_comp,η0,η,dη_sum)
             end
 
         elseif time_scheme == "SSPRK2"  # s-stage 2nd order SSPRK
@@ -116,6 +145,12 @@ function time_integration(  Prog::PrognosticVars{Tprog},
 
             @unpack s,kn,mn,kna,knb,Δt_Δnc,Δt_Δn = S.constants.SSPRK3c
 
+            # if compensated
+            #     fill!(du_sum,zero(Tprog))
+            #     fill!(dv_sum,zero(Tprog))
+            #     fill!(dη_sum,zero(Tprog))
+            # end
+
             for rki = 2:s+1       # number of stages (from 2:s+1 to match Ketcheson et al 2014)
                 if rki > 2
                     ghost_points_η!(η1,S)
@@ -134,6 +169,11 @@ function time_integration(  Prog::PrognosticVars{Tprog},
                 else                                # normal update case
                     axb!(u1,Δt_Δn,du)   
                     axb!(v1,Δt_Δn,dv)
+
+                    # if compensated
+                    #     axb!(du_sum,Δt_Δn,du)   
+                    #     axb!(dv_sum,Δt_Δn,dv)
+                    # end
                 end
 
                 # semi-implicit for continuity equation, use new u1,v1 to calcualte dη
@@ -146,6 +186,9 @@ function time_integration(  Prog::PrognosticVars{Tprog},
                     dxaybzc!(η1,kna,η1,knb,η0,Δt_Δnc,dη)
                 else
                     axb!(η1,Δt_Δn,dη)
+                    # if compensated
+                    #     axb!(dη_sum,Δt_Δn,dη)
+                    # end
                 end
 
                 # special stage that is needed later for the kn-th stage, store in u0,v0,η0 therefore
@@ -253,14 +296,15 @@ function time_integration(  Prog::PrognosticVars{Tprog},
 end
 
 """Add to a x multiplied with b. a += x*b """
-function axb!(a::Array{T,2},x::T,b::Array{T,2}) where {T<:AbstractFloat}
+function axb!(a::Matrix{T},x::Real,b::Matrix{T}) where {T<:AbstractFloat}
     m,n = size(a)
     @boundscheck (m,n) == size(b) || throw(BoundsError())
 
-    #TODO @simd?
+    xT = convert(T,x)
+
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-           a[i,j] += x*b[i,j]
+           a[i,j] += xT*b[i,j]
         end
     end
 end
@@ -278,13 +322,28 @@ function caxb!(c::Array{T,2},a::Array{T,2},x::T,b::Array{T,2}) where {T<:Abstrac
     end
 end
 
+"""d equals add a minus b minus c. c = (a - b) - c."""
+function dambmc!(d::Matrix{T},a::Matrix{T},b::Matrix{T},c::Matrix{T}) where {T<:AbstractFloat}
+    m,n = size(a)
+    @boundscheck (m,n) == size(b) || throw(BoundsError())
+    @boundscheck (m,n) == size(c) || throw(BoundsError())
+    @boundscheck (m,n) == size(d) || throw(BoundsError())
+
+    @inbounds for j ∈ 1:n
+        for i ∈ 1:m
+           d[i,j] = (a[i,j] - b[i,j]) - c[i,j]
+        end
+    end
+end
+
 """c equals add x multiplied to a plus b. c = x*(a+b) """
 function cxab!(c::Array{T,2},x::Real,a::Array{T,2},b::Array{T,2}) where {T<:AbstractFloat}
     m,n = size(a)
     @boundscheck (m,n) == size(b) || throw(BoundsError())
     @boundscheck (m,n) == size(c) || throw(BoundsError())
 
-    xT = T(x)
+    xT = convert(T,x)
+
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
            c[i,j] = xT*(a[i,j] + b[i,j])
@@ -298,8 +357,9 @@ function cxayb!(c::Array{T,2},x::Real,a::Array{T,2},y::Real,b::Array{T,2}) where
     @boundscheck (m,n) == size(b) || throw(BoundsError())
     @boundscheck (m,n) == size(c) || throw(BoundsError())
 
-    xT = T(x)
-    yT = T(y)
+    xT = convert(T,x)
+    yT = convert(T,y)
+
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
            c[i,j] = xT*a[i,j] + yT*b[i,j]
@@ -317,9 +377,9 @@ function dxaybzc!(  d::Array{T,2},
     @boundscheck (m,n) == size(c) || throw(BoundsError())
     @boundscheck (m,n) == size(d) || throw(BoundsError())
 
-    xT = T(x)       # convert to type T
-    yT = T(y)
-    zT = T(z)
+    xT = convert(T,x)
+    yT = convert(T,y)
+    zT = convert(T,z)
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m