Merge pull request #113 from milankl/mixedprec

Mixedprec

Author: Milan K

.appveyor.yml

@@ -1,7 +1,7 @@
 environment:
   matrix:
   - julia_version: 1
-  - julia_version: 1.2
+  - julia_version: 1.3
   - julia_version: nightly
 
 platform:

.cirrus.yml

@@ -5,7 +5,7 @@ task:
   env:
     matrix:
       - JULIA_VERSION: 1.0
-      - JULIA_VERSION: 1.2
+      - JULIA_VERSION: 1.3
       - JULIA_VERSION: nightly
   install_script:
     - sh -c "$(fetch https://raw.githubusercontent.com/ararslan/CirrusCI.jl/master/bin/install.sh -o -)"

.travis.yml

@@ -5,7 +5,7 @@ os:
   - osx
 julia:
   - 1.0
-  - 1.2
+  - 1.3
   - nightly
 notifications:
   email: false

src/Constants.jl

@@ -1,11 +1,11 @@
-struct Constants{T<:AbstractFloat}
+struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat}
 
     # RUNGE-KUTTA COEFFICIENTS 3rd/4th order including timestep Δt
-    RKaΔt::Array{T,1}
-    RKbΔt::Array{T,1}
+    RKaΔt::Array{Tprog,1}
+    RKbΔt::Array{Tprog,1}
 
     # BOUNDARY CONDITIONS
-    one_minus_α::T      # tangential boundary condition for the ghost-point copy
+    one_minus_α::Tprog      # tangential boundary condition for the ghost-point copy
 
     # PHYSICAL CONSTANTS
     g::T                    # gravity
@@ -17,23 +17,23 @@ struct Constants{T<:AbstractFloat}
     rSST::T                 # tracer restoring timescale
     jSST::T                 # tracer consumption timescale
     SSTmin::T               # tracer minimum
-    ωyr::Float64            # # frequency [1/s] of Kelvin pumping (including 2π)
+    ωyr::Float64            # frequency [1/s] of Kelvin pumping (including 2π)
 end
 
 """Generator function for the mutable struct Constants."""
-function Constants{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
+function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<:AbstractFloat}
 
     # Runge-Kutta 3rd/4th order coefficients including time step Δt
     # (which includes the grid spacing Δ too)
     if P.RKo == 3     # version 2
-        RKaΔt = T.([1/4,0.,3/4]*G.dtint/G.Δ)
-        RKbΔt = T.([1/3,2/3]*G.dtint/G.Δ)
+        RKaΔt = Tprog.([1/4,0.,3/4]*G.dtint/G.Δ)
+        RKbΔt = Tprog.([1/3,2/3]*G.dtint/G.Δ)
     elseif P.RKo == 4
-        RKaΔt = T.([1/6,1/3,1/3,1/6]*G.dtint/G.Δ)
-        RKbΔt = T.([.5,.5,1.]*G.dtint/G.Δ)
+        RKaΔt = Tprog.([1/6,1/3,1/3,1/6]*G.dtint/G.Δ)
+        RKbΔt = Tprog.([.5,.5,1.]*G.dtint/G.Δ)
     end
 
-    one_minus_α = T(1-P.α)    # for the ghost point copy/tangential boundary conditions
+    one_minus_α = Tprog(1-P.α)    # for the ghost point copy/tangential boundary conditions
     g = T(P.g)                # gravity - for Bernoulli potential
 
     # BOTTOM FRICTION COEFFICENTS
@@ -57,5 +57,5 @@ function Constants{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
     # SURFACE FORCING
     ωyr = -2π*P.ωyr/24/365/3600
 
-    return Constants{T}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin,ωyr)
+    return Constants{T,Tprog}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin,ωyr)
 end

src/Continuity.jl

@@ -1,9 +1,9 @@
 """Continuity equation's right-hand side with surface relaxation
 -∂x(uh) - ∂y(vh) + γ*(η_ref - η)."""
-function continuity_surf_relax!(η::AbstractMatrix,
-                                Diag::DiagnosticVars,
-                                S::ModelSetup,
-                                t::Int) where {T<:AbstractFloat}
+function continuity_surf_relax!(η::Array{T,2},
+                                Diag::DiagnosticVars{T,Tprog},
+                                S::ModelSetup{T,Tprog},
+                                t::Int) where {T,Tprog}
 
     @unpack dη = Diag.Tendencies
     @unpack dUdx,dVdy = Diag.VolumeFluxes
@@ -18,17 +18,15 @@ function continuity_surf_relax!(η::AbstractMatrix,
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            dη[i+1,j+1] = -(dUdx[i,j+1] + dVdy[i+1,j]) + γ*(η_ref[i,j]-η[i+1,j+1])
+            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(γ*(η_ref[i,j]-η[i+1,j+1]))
         end
     end
 end
 
 """Continuity equation's right-hand side with time&space dependent forcing."""
-function continuity_forcing!(   ::Type{T},
-                                η::AbstractMatrix,
-                                Diag::DiagnosticVars,
-                                S::ModelSetup,
-                                t::Int) where {T<:AbstractFloat}
+function continuity_forcing!(   Diag::DiagnosticVars{T,Tprog},
+                                S::ModelSetup{T,Tprog},
+                                t::Int) where {T,Tprog}
 
     @unpack dη = Diag.Tendencies
     @unpack dUdx,dVdy = Diag.VolumeFluxes
@@ -45,17 +43,15 @@ function continuity_forcing!(   ::Type{T},
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            dη[i+1,j+1] = -(dUdx[i,j+1] + dVdy[i+1,j]) + Fηtt*Fη[i,j]
+            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(Fηtt*Fη[i,j])
         end
     end
 end
 
 """Continuity equation's right-hand side -∂x(uh) - ∂y(vh) without forcing."""
-function continuity_itself!(::Type{T},
-                            η::AbstractMatrix,
-                            Diag::DiagnosticVars,
-                            S::ModelSetup,
-                            t::Int) where {T<:AbstractFloat}
+function continuity_itself!(Diag::DiagnosticVars{T,Tprog},
+                            S::ModelSetup{T,Tprog},
+                            t::Int) where {T,Tprog}
 
     @unpack dη = Diag.Tendencies
     @unpack dUdx,dVdy = Diag.VolumeFluxes
@@ -66,23 +62,23 @@ function continuity_itself!(::Type{T},
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            dη[i+1,j+1] = -(dUdx[i,j+1] + dVdy[i+1,j])
+            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j]))
         end
     end
 end
 
 
 """Transit function to call the specified continuity function."""
-function continuity!(   η::Array{T,2},
+function continuity!(   η::AbstractMatrix,
                         Diag::DiagnosticVars,
                         S::ModelSetup,
-                        t::Int) where {T<:AbstractFloat}
+                        t::Int)
 
     if S.parameters.surface_relax
-        continuity_surf_relax!(T,η,Diag,S,t)
+        continuity_surf_relax!(η,Diag,S,t)
     elseif S.parameters.surface_forcing
-        continuity_forcing!(T,η,Diag,S,t)
+        continuity_forcing!(Diag,S,t)
     else
-        continuity_itself!(T,η,Diag,S,t)
+        continuity_itself!(Diag,S,t)
     end
 end

src/DefaultParameters.jl

@@ -2,6 +2,9 @@
 
     T::DataType=Float32                 # number format
 
+    Tprog::DataType=T                   # number format for prognostic variables
+    Tcomm::DataType=T                   # number format for ghost-point copies
+
     # DOMAIN RESOLUTION AND RATIO
     nx::Int=100                         # number of grid cells in x-direction
     Lx::Real=2000e3                     # length of the domain in x-direction [m]
@@ -30,7 +33,7 @@
     surface_relax::Bool=false           # yes?
     t_relax::Real=100.                  # time scale of the relaxation [days]
     η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
-    η_refw::Real=50e3                  # width [m] of the tangent used for the interface relaxation
+    η_refw::Real=50e3                   # width [m] of the tangent used for the interface relaxation
 
     # SURFACE FORCING (Currently only Kelvin wave pumping at Eq.)
     surface_forcing::Bool=false         # yes?
@@ -73,11 +76,11 @@
     sst_initial::String="south"         # "west", "south", "rect", "flat" or "restart"
     sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                         # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
-    Uadv::Real=0.15                     # Velocity scale [m/s] for tracer advection
+    Uadv::Real=0.25                     # Velocity scale [m/s] for tracer advection
     SSTmax::Real=1.                     # tracer (sea surface temperature) max for restoring
     SSTmin::Real=0.                     # tracer (sea surface temperature) min for restoring
     τSST::Real=500.                     # tracer restoring time scale [days]
-    jSST::Real=50*365.                  # tracer consumption [days]
+    jSST::Real=365.                     # tracer consumption [days]
     SST_λ0::Real=222e3                  # [m] transition position of relaxation timescale
     SST_λs::Real=111e3                  # [m] transition width of relaxation timescale
     SST_γ0::Real=8.35                   # [days] injection time scale

src/Diffusion.jl

@@ -218,10 +218,10 @@ function viscous_tensor_constant!(  Diag::DiagnosticVars,
 end
 
 """Update u with bottom friction tendency (Bu,Bv) and biharmonic viscosity."""
-function add_drag_diff_tendencies!( u::AbstractMatrix,
-                                    v::AbstractMatrix,
-                                    Diag::DiagnosticVars,
-                                    S::ModelSetup)
+function add_drag_diff_tendencies!( u::Array{Tprog,2},
+                                    v::Array{Tprog,2},
+                                    Diag::DiagnosticVars{T,Tprog},
+                                    S::ModelSetup{T,Tprog}) where {T,Tprog}
 
     @unpack Bu,Bv = Diag.Bottomdrag
     @unpack LLu1,LLu2,LLv1,LLv2 = Diag.Smagorinsky
@@ -234,7 +234,7 @@ function add_drag_diff_tendencies!( u::AbstractMatrix,
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            u[i+2,j+2] += Δt_diff*(Bu[i+1-ep,j+1] + LLu1[i,j+1] + LLu2[i+1-ep,j])
+            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
 
@@ -245,7 +245,7 @@ function add_drag_diff_tendencies!( u::AbstractMatrix,
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-             v[i+2,j+2] += Δt_diff*(Bv[i+1,j+1] + LLv1[i,j+1] + LLv2[i+1,j])
+             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

src/Grid.jl

@@ -1,4 +1,4 @@
-@with_kw struct Grid{T<:AbstractFloat}
+@with_kw struct Grid{T<:AbstractFloat,Tprog<:AbstractFloat}
 
     # Parameters taken from Parameter struct
     nx::Int                             # number of grid cells in x-direction
@@ -71,7 +71,7 @@
     nt::Int = Int(ceil(Ndays*3600*24/dtint))    # number of time steps to integrate
     dt::T = T(dtint)                            # time step [s]
     Δt::T = T(dtint/Δ)                          # time step divided by grid spacing [s/m]
-    Δt_diff::T = T(nstep_diff*dtint/Δ)          # time step for diffusive terms
+    Δt_diff::Tprog = Tprog(nstep_diff*dtint/Δ)  # time step for diffusive terms
 
     # TIME STEPS FOR ADVECTION
     nadvstep::Int = max(1,Int(floor(Δ/Uadv/dtint)))         # advection each n time steps
@@ -80,7 +80,7 @@
     dtadvu::T = T(dtadvint*nx/Lx)                           # Rescaled advection time step for u [s/m]
     dtadvv::T = T(dtadvint*ny/Ly)                           # Rescaled advection time step for v [s/m]
     half_dtadvu::T = T(dtadvint*nx/Lx/2)                    # dtadvu/2
-    half_dtadvv::T = T(dtadvint*ny/Ly/2)                    # dtadvv/2 
+    half_dtadvv::T = T(dtadvint*ny/Ly/2)                    # dtadvv/2
 
     # N TIME STEPS FOR OUTPUT
     nout::Int = Int(floor(output_dt*3600/dtint))            # output every n time steps
@@ -99,24 +99,24 @@ end
 function yy_q(bc::String,x_q_halo::AbstractVector,y_q_halo::AbstractVector)
     if bc == "periodic"
         # points on the right edge needed too
-        xx_q,yy_q = meshgrid(x_q_halo[3:end-1],y_q_halo[3:end-2])
+        _,yy_q = meshgrid(x_q_halo[3:end-1],y_q_halo[3:end-2])
     else
-        xx_q,yy_q = meshgrid(x_q_halo[3:end-2],y_q_halo[3:end-2])
+        _,yy_q = meshgrid(x_q_halo[3:end-2],y_q_halo[3:end-2])
     end
 
     return yy_q
 end
 
 
 """Generator function for the Grid struct."""
-function Grid{T}(P::Parameter) where {T<:AbstractFloat}
+function Grid{T,Tprog}(P::Parameter) where {T<:AbstractFloat,Tprog<:AbstractFloat}
     @unpack nx,Lx,L_ratio = P
     @unpack bc,g,H,cfl = P
     @unpack Ndays,nstep_diff,nstep_advcor = P
     @unpack Uadv,output_dt = P
     @unpack ϕ,ω,R = P
 
-    return Grid{T}(nx=nx,Lx=Lx,L_ratio=L_ratio,bc=bc,g=g,H=H,cfl=cfl,Ndays=Ndays,
+    return Grid{T,Tprog}(nx=nx,Lx=Lx,L_ratio=L_ratio,bc=bc,g=g,H=H,cfl=cfl,Ndays=Ndays,
                 nstep_diff=nstep_diff,nstep_advcor=nstep_advcor,Uadv=Uadv,output_dt=output_dt,
                 ϕ=ϕ,ω=ω,R=R)
 end

src/InitialConditions.jl

@@ -10,6 +10,19 @@ struct PrognosticVars{T<:AbstractFloat}
     sst::Array{T,2}         # tracer / sea surface temperature
 end
 
+"""Zero generator function for Grid G as argument."""
+function PrognosticVars{T}(G::Grid) where {T<:AbstractFloat}
+    @unpack nux,nuy,nvx,nvy,nx,ny = G
+    @unpack halo,haloη,halosstx,halossty = G
+
+    u = zeros(T,nux+2halo,nuy+2halo)
+    v = zeros(T,nvx+2halo,nvy+2halo)
+    η = zeros(T,nx+2haloη,ny+2haloη)
+    sst = zeros(T,nx+2halosstx,ny+2halossty)
+
+    return PrognosticVars{T}(u,v,η,sst)
+end
+
 function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
 
     ## PROGNOSTIC VARIABLES U,V,η
@@ -49,7 +62,7 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
         NetCDF.close(ncη)
 
         # remove singleton time dimension
-        # and convert from Float32 to Numtype
+        # and convert from Float32 to type T
         u = T.(reshape(u,size(u)[1:2]))
         v = T.(reshape(v,size(v)[1:2]))
         η = T.(reshape(η,size(η)[1:2]))

src/ModelSetup.jl

@@ -1,6 +1,6 @@
-struct ModelSetup{T<:AbstractFloat}
+struct ModelSetup{T<:AbstractFloat,Tprog<:AbstractFloat}
     parameters::Parameter
-    grid::Grid{T}
-    constants::Constants{T}
+    grid::Grid{T,Tprog}
+    constants::Constants{T,Tprog}
     forcing::Forcing{T}
 end