Merge pull request #155 from milankl/scale

scale also sst

Author: milankl

README.md

@@ -1,7 +1,7 @@
 # ShallowWaters.jl - A type-flexible 16-bit shallow water model
 [![CI](https://github.com/milankl/ShallowWaters.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/milankl/ShallowWaters.jl/actions/workflows/CI.yml)
 [![DOI](https://zenodo.org/badge/132787050.svg)](https://zenodo.org/badge/latestdoi/132787050)
-![sst](figs/sst_posit16.png?raw=true "SST")
+![sst](figs/isambard_float16.png?raw=true "Float16 simulation with ShallowWaters.jl on Isambard's A64FX")
 
 A shallow water model with a focus on type-flexibility and 16-bit number formats. ShallowWaters allows for Float64/32/16, 
 [Posit32/16/8](https://github.com/milankl/SoftPosit.jl), [BFloat16](https://github.com/JuliaComputing/BFloat16s.jl), 

figs/isambard_float16.png

@@ -4,6 +4,7 @@
 
     Tprog=T                   # number format for prognostic variables
     Tcomm=Tprog               # number format for ghost-point copies
+    Tini=Tprog                # number format to reduce precision for initial conditions
 
     # DOMAIN RESOLUTION AND RATIO
     nx::Int=100                         # number of grid cells in x-direction
@@ -157,6 +158,7 @@ Creates a Parameter struct with following options and default values
 
     Tprog=T                   # number format for prognostic variables
     Tcomm=Tprog               # number format for ghost-point copies
+    Tini=Tprog                # number format to reduce precision for initial conditions
 
     # DOMAIN RESOLUTION AND RATIO
     nx::Int=100                         # number of grid cells in x-direction

src/default_parameters.jl

@@ -173,4 +173,4 @@ function progress!(feedback::Feedback)
             flush(progress_txt)
         end
     end
-end
+end

src/feedback.jl

@@ -25,6 +25,7 @@ function add_halo(  u::Array{T,2},
     @unpack scale = S.constants
     u *= scale
     v *= scale
+    sst *= scale
 
     ghost_points!(u,v,η,S)
     ghost_points_sst!(sst,S)
@@ -41,10 +42,11 @@ function remove_halo(   u::Array{T,2},
     @unpack halo,haloη,halosstx,halossty = S.grid
     @unpack scale_inv = S.constants
 
+    # undo scaling as well
     ucut = scale_inv*u[halo+1:end-halo,halo+1:end-halo]
     vcut = scale_inv*v[halo+1:end-halo,halo+1:end-halo]
     ηcut = η[haloη+1:end-haloη,haloη+1:end-haloη]
-    sstcut = sst[halosstx+1:end-halosstx,halossty+1:end-halossty]
+    sstcut = scale_inv*sst[halosstx+1:end-halosstx,halossty+1:end-halossty]
 
     return ucut,vcut,ηcut,sstcut
 end

src/ghost_points.jl

@@ -1,13 +1,11 @@
 """Calculates the 2nd order centred gradient in x-direction on any grid (u,v,T or q).
 The size of dudx must be m-1,n compared to m,n = size(u)"""
-function ∂x!(∂ₓu::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
-    m,n = size(∂ₓu)
+function ∂x!(dudx::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
+    m,n = size(dudx)
     @boundscheck (m+1,n) == size(u) || throw(BoundsError())
 
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            ∂ₓu[i,j] = u[i+1,j] - u[i,j]
-        end
+    @inbounds for j ∈ 1:n, i ∈ 1:m
+        dudx[i,j] = u[i+1,j] - u[i,j]
     end
 end
 
@@ -17,16 +15,14 @@ function ∂y!(dudy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
     m,n = size(dudy)
     @boundscheck (m,n+1) == size(u) || throw(BoundsError())
 
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
+    @inbounds for j ∈ 1:n, i ∈ 1:m
             dudy[i,j] = u[i,j+1] - u[i,j]
-        end
     end
 end
 
 """ ∇² is the 2nd order centred Laplace-operator ∂/∂x^2 + ∂/∂y^2.
 The 1/Δ²-factor is omitted and moved into the viscosity coefficient."""
-function ∇²!(du::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
+function ∇²!(du::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
     m, n = size(du)
     @boundscheck (m+2,n+2) == size(u) || throw(BoundsError())
 
@@ -91,4 +87,4 @@ function ∇²(u::Array{T,2},Δ::Real=1) where {T<:AbstractFloat}
         end
     end
     return du
-end
+end

src/gradients.jl

@@ -26,9 +26,9 @@ end
 function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
 
     ## PROGNOSTIC VARIABLES U,V,η
-
     @unpack nux,nuy,nvx,nvy,nx,ny = S.grid
     @unpack initial_cond = S.parameters
+    @unpack Tini = S.parameters
 
     if initial_cond == "rest"
 
@@ -53,15 +53,12 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
         end
 
         u = ncu.vars["u"][:,:,init_starti]
-        # NetCDF.close(ncu)
 
         ncv = NetCDF.open(joinpath(inirunpath,"v.nc"))
         v = ncv.vars["v"][:,:,init_starti]
-        # NetCDF.close(ncv)
 
         ncη = NetCDF.open(joinpath(inirunpath,"eta.nc"))
         η = ncη.vars["eta"][:,:,init_starti]
-        # NetCDF.close(ncη)
 
         # remove singleton time dimension
         u = reshape(u,size(u)[1:2])
@@ -122,7 +119,7 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
     ## SST
 
     @unpack SSTmin, SSTmax, SSTw, SSTϕ = S.parameters
-    @unpack sst_initial = S.parameters
+    @unpack sst_initial,scale = S.parameters
     @unpack x_T,y_T,Lx,Ly = S.grid
 
     xx_T,yy_T = meshgrid(x_T,y_T)
@@ -145,19 +142,18 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
     if initial_cond == "ncfile" && sst_initial == "restart"
         ncsst = NetCDF.open(joinpath(inirunpath,"sst.nc"))
         sst = ncsst.vars["sst"][:,:,init_starti]
-        # NetCDF.close(ncsst)
 
         sst = reshape(sst,size(sst)[1:2])
     end
 
     # Convert to number format T
-    sst = T.(sst)
-    u = T.(u)
-    v = T.(v)
-    η = T.(η)
+    # allow for comparable initial conditions via Tini
+    sst = T.(Tini.(sst))
+    u = T.(Tini.(u))
+    v = T.(Tini.(v))
+    η = T.(Tini.(η))
 
     #TODO SST INTERPOLATION
-
     u,v,η,sst = add_halo(u,v,η,sst,S)
 
     return PrognosticVars{T}(u,v,η,sst)

src/initial_conditions.jl

@@ -1,6 +1,6 @@
 """Linear interpolation of a variable u in the x-direction.
 m,n = size(ux) must be m+1,n = size(u)."""
-function Ix!(ux::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
+function Ix!(ux::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
     m, n = size(ux)
     @boundscheck (m+1,n) == size(u) || throw(BoundsError())
 
@@ -15,7 +15,7 @@ end
 
 """ Linear interpolation a variable u in the y-direction.
     m,n = size(uy) must be m,n+1 = size(u)."""
-function Iy!(uy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
+function Iy!(uy::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
     m,n = size(uy)
     @boundscheck (m,n+1) == size(u) || throw(BoundsError())
 
@@ -36,11 +36,9 @@ function Ixy!(uxy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
 
     one_quarter = convert(T,0.25)
 
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            uxy[i,j] = one_quarter*(u[i,j] + u[i+1,j]) + 
-                    one_quarter*(u[i,j+1] + u[i+1,j+1])
-        end
+    @inbounds for j in 1:n, i in 1:m
+        uxy[i,j] = one_quarter*(u[i,j] + u[i+1,j]) + 
+                one_quarter*(u[i,j+1] + u[i+1,j+1])
     end
 end
 

src/interpolations.jl

@@ -105,7 +105,7 @@ function output_nc!(i::Int,
             NetCDF.putvar(ncs.η,"eta",η,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.sst != nothing
-            @views sst = Float32.(Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
+            @views sst = Float32.(scale_inv*Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
             NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.q != nothing