Skip to content

Commit f163176

Browse files
milanklmilankl
committed
scale also sst
1 parent 03faffb commit f163176

File tree

6 files changed

+75
-26
lines changed

6 files changed

+75
-26
lines changed

‎src/feedback.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -173,4 +173,4 @@ function progress!(feedback::Feedback)
173173
flush(progress_txt)
174174
end
175175
end
176-
end
176+
end

‎src/ghost_points.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ function add_halo( u::Array{T,2},
2525
@unpack scale = S.constants
2626
u *= scale
2727
v *= scale
28+
sst *= scale
2829

2930
ghost_points!(u,v,η,S)
3031
ghost_points_sst!(sst,S)

‎src/gradients.jl

Lines changed: 66 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,11 @@
11
"""Calculates the 2nd order centred gradient in x-direction on any grid (u,v,T or q).
22
The size of dudx must be m-1,n compared to m,n = size(u)"""
3-
function ∂x!(∂ₓu::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
4-
m,n = size(∂ₓu)
3+
function ∂x!(dudx::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
4+
m,n = size(dudx)
55
@boundscheck (m+1,n) == size(u) || throw(BoundsError())
66

7-
@inbounds for j ∈ 1:n
8-
for i ∈ 1:m
9-
∂ₓu[i,j] = u[i+1,j] - u[i,j]
10-
end
7+
@inbounds for j ∈ 1:n, i ∈ 1:m
8+
dudx[i,j] = u[i+1,j] - u[i,j]
119
end
1210
end
1311

@@ -17,16 +15,14 @@ function ∂y!(dudy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
1715
m,n = size(dudy)
1816
@boundscheck (m,n+1) == size(u) || throw(BoundsError())
1917

20-
@inbounds for j ∈ 1:n
21-
for i ∈ 1:m
18+
@inbounds for j ∈ 1:n, i ∈ 1:m
2219
dudy[i,j] = u[i,j+1] - u[i,j]
23-
end
2420
end
2521
end
2622

2723
""" ∇² is the 2nd order centred Laplace-operator ∂/∂x^2 + ∂/∂y^2.
2824
The 1/Δ²-factor is omitted and moved into the viscosity coefficient."""
29-
function ∇²!(du::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
25+
function ∇²!(du::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
3026
m, n = size(du)
3127
@boundscheck (m+2,n+2) == size(u) || throw(BoundsError())
3228

@@ -92,3 +88,63 @@ function ∇²(u::Array{T,2},Δ::Real=1) where {T<:AbstractFloat}
9288
end
9389
return du
9490
end
91+
92+
93+
function benchmark_it(N::Vector{Int})
94+
95+
timings = fill(0.0,3,length(N))
96+
97+
for (i,n) in enumerate(N)
98+
for (j,T) in enumerate([Float64,Float32,Float16])
99+
A = rand(T,n,n)
100+
B = rand(T,n,n)
101+
C = rand(T,n,n)
102+
D = rand(T,n,n)
103+
E = rand(T,n,n)
104+
F = rand(T,n,n)
105+
R = zero(A)
106+
timings[j,i] = minimum(@benchmark add!($R,$A,$B,$C,$D,$E,$F) samples=100 evals=1).time
107+
end
108+
end
109+
110+
print("N= ")
111+
for n in N
112+
print(@sprintf(" %4d,",n))
113+
end
114+
println()
115+
116+
for (iT,T) in enumerate([Float64,Float32,Float16])
117+
print("$T:")
118+
for t in 1:length(N)
119+
trel = timings[1,t]/timings[iT,t]
120+
print(@sprintf(" %.3fx,",trel))
121+
end
122+
println()
123+
end
124+
end
125+
126+
127+
function add!(R::Matrix{T},
128+
A::Matrix{T},
129+
B::Matrix{T},
130+
C::Matrix{T},
131+
D::Matrix{T},
132+
E::Matrix{T},
133+
F::Matrix{T},
134+
G::Matrix{T},
135+
H::Matrix{T}) where {T<:AbstractFloat}
136+
137+
m,n = size(R)
138+
@boundscheck size(R) == size(A) || throw(BoundsError())
139+
@boundscheck size(R) == size(B) || throw(BoundsError())
140+
@boundscheck size(R) == size(C) || throw(BoundsError())
141+
@boundscheck size(R) == size(D) || throw(BoundsError())
142+
@boundscheck size(R) == size(E) || throw(BoundsError())
143+
@boundscheck size(R) == size(F) || throw(BoundsError())
144+
@boundscheck size(R) == size(G) || throw(BoundsError())
145+
@boundscheck size(R) == size(H) || throw(BoundsError())
146+
147+
@inbounds for i in eachindex(R)
148+
R[i] = A[i] + B[i] + C[i] + D[i] + E[i] + F[i] + G[i] + H[i]
149+
end
150+
end

‎src/initial_conditions.jl

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@ end
2626
function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
2727

2828
## PROGNOSTIC VARIABLES U,V,η
29-
3029
@unpack nux,nuy,nvx,nvy,nx,ny = S.grid
3130
@unpack initial_cond = S.parameters
3231

@@ -53,15 +52,12 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
5352
end
5453

5554
u = ncu.vars["u"][:,:,init_starti]
56-
# NetCDF.close(ncu)
5755

5856
ncv = NetCDF.open(joinpath(inirunpath,"v.nc"))
5957
v = ncv.vars["v"][:,:,init_starti]
60-
# NetCDF.close(ncv)
6158

6259
ncη = NetCDF.open(joinpath(inirunpath,"eta.nc"))
6360
η = ncη.vars["eta"][:,:,init_starti]
64-
# NetCDF.close(ncη)
6561

6662
# remove singleton time dimension
6763
u = reshape(u,size(u)[1:2])
@@ -122,7 +118,7 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
122118
## SST
123119

124120
@unpack SSTmin, SSTmax, SSTw, SSTϕ = S.parameters
125-
@unpack sst_initial = S.parameters
121+
@unpack sst_initial,scale = S.parameters
126122
@unpack x_T,y_T,Lx,Ly = S.grid
127123

128124
xx_T,yy_T = meshgrid(x_T,y_T)
@@ -145,7 +141,6 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
145141
if initial_cond == "ncfile" && sst_initial == "restart"
146142
ncsst = NetCDF.open(joinpath(inirunpath,"sst.nc"))
147143
sst = ncsst.vars["sst"][:,:,init_starti]
148-
# NetCDF.close(ncsst)
149144

150145
sst = reshape(sst,size(sst)[1:2])
151146
end
@@ -157,7 +152,6 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
157152
η = T.(η)
158153

159154
#TODO SST INTERPOLATION
160-
161155
u,v,η,sst = add_halo(u,v,η,sst,S)
162156

163157
return PrognosticVars{T}(u,v,η,sst)

‎src/interpolations.jl

Lines changed: 5 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
"""Linear interpolation of a variable u in the x-direction.
22
m,n = size(ux) must be m+1,n = size(u)."""
3-
function Ix!(ux::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
3+
function Ix!(ux::Matrix{T},u::Matrix{T}) where {T<:AbstractFloat}
44
m, n = size(ux)
55
@boundscheck (m+1,n) == size(u) || throw(BoundsError())
66

@@ -15,7 +15,7 @@ end
1515

1616
""" Linear interpolation a variable u in the y-direction.
1717
m,n = size(uy) must be m,n+1 = size(u)."""
18-
function Iy!(uy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
18+
function Iy!(uy::Matrix{T},u::Matrix{2}) where {T<:AbstractFloat}
1919
m,n = size(uy)
2020
@boundscheck (m,n+1) == size(u) || throw(BoundsError())
2121

@@ -36,11 +36,9 @@ function Ixy!(uxy::Array{T,2},u::Array{T,2}) where {T<:AbstractFloat}
3636

3737
one_quarter = convert(T,0.25)
3838

39-
@inbounds for j ∈ 1:n
40-
for i ∈ 1:m
41-
uxy[i,j] = one_quarter*(u[i,j] + u[i+1,j]) +
42-
one_quarter*(u[i,j+1] + u[i+1,j+1])
43-
end
39+
@inbounds for j in 1:n, i in 1:m
40+
uxy[i,j] = one_quarter*(u[i,j] + u[i+1,j]) +
41+
one_quarter*(u[i,j+1] + u[i+1,j+1])
4442
end
4543
end
4644

‎src/output.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,7 @@ function output_nc!(i::Int,
105105
NetCDF.putvar(ncs.η,"eta",η,start=[1,1,iout],count=[-1,-1,1])
106106
end
107107
if ncs.sst != nothing
108-
@views sst = Float32.(Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
108+
@views sst = Float32.(scale_inv*Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
109109
NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
110110
end
111111
if ncs.q != nothing

0 commit comments

Comments
 (0)