Merge pull request #144 from milankl/rk2

Rk2

Author: milankl

.appveyor.yml

@@ -1,6 +1,7 @@
 environment:
   matrix:
-  - julia_version: 1.4
+  - julia_version: 1.2
+  - julia_version: 1.5
   - julia_version: nightly
 
 platform:

.travis.yml

@@ -4,7 +4,8 @@ os:
   - linux
   - osx
 julia:
-  - 1.4
+  - 1.2
+  - 1.5
   - nightly
 notifications:
   email: false

Project.toml

@@ -1,7 +1,7 @@
 name = "ShallowWaters"
 uuid = "56019723-2d87-4a65-81ff-59d5d8913e3c"
 authors = ["Milan Kloewer"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -14,7 +14,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Interpolations = "^0.12"
 NetCDF = "^0.10"
 Parameters = "^0.12"
-julia = "1"
+julia = "^1.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Constants.jl

@@ -25,9 +25,13 @@ end
 """Generator function for the mutable struct Constants."""
 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
+    # Runge-Kutta 2nd/3rd/4th order coefficients including time step Δt and grid spacing Δ
+    # a are the coefficents to sum the rhs on the fly, such that sum=1
+    # b are the coefficents for the update that used for a new evaluation of the RHS
+    if P.RKo == 2     # Heun's method
+        RKaΔt = Tprog.([1/2,1/2]*G.dtint/G.Δ)
+        RKbΔt = Tprog.([1]*G.dtint/G.Δ)
+    elseif P.RKo == 3     # version 2 / Heun's 3rd order
         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

src/DefaultParameters.jl

@@ -95,7 +95,7 @@
 
     # OUTPUT OPTIONS
     output::Bool=false                  # netcdf output?
-    output_vars::Array{String,1}=["u","v","η","sst","q","ζ"]  # which variables to output? "du","dv","dη" also allowed
+    output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? "du","dv","dη" also allowed
     output_dt::Real=6                   # output time step [hours]
     outpath::String=pwd()               # path to output folder
 
@@ -118,7 +118,7 @@
     @assert topo_width > 0.0    "topo_width has to be >0, $topo_width given."
     @assert t_relax > 0.0       "t_relax has to be >0, $t_relax given."
     @assert η_refw > 0.0        "η_refw has to be >0, $η_refw given."
-    @assert RKo in [3,4]        "RKo has to be 3 or 4, $RKo given."
+    @assert RKo in [2,3,4]        "RKo has to be 2,3 or 4, $RKo given."
     @assert Ndays > 0.0         "Ndays has to be >0, $Ndays given."
     @assert nstep_diff > 0      "nstep_diff has to be >0, $nstep_diff given."
     @assert nstep_advcor >= 0   "nstep_advcor has to be >=0, $nstep_advcor given."
@@ -191,8 +191,8 @@ Creates a Parameter struct with following options and default values
     wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing
 
     # TIME STEPPING OPTIONS
-    RKo::Int=4                          # Order of the RK time stepping scheme (3 or 4)
-    cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
+    RKo::Int=4                          # Order of the RK time stepping scheme (2,3 or 4)
+    cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3, 0.1 for RK2)
     Ndays::Real=10.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.
@@ -239,7 +239,7 @@ Creates a Parameter struct with following options and default values
 
     # OUTPUT OPTIONS
     output::Bool=false                  # netcdf output?
-    output_vars::Array{String,1}=["u","v","η","sst","q","ζ"]  # which variables to output? "du","dv","dη" also allowed
+    output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? q,ζ,du,dv,dη also allowed.
     output_dt::Real=6                   # output time step [hours]
     outpath::String=pwd()               # path to output folder
 

src/Output.jl

@@ -65,7 +65,8 @@ function nc_create( x::Array{T,1},
     tvar = NcVar("t",tdim,t=Int32)
 
     nc = NetCDF.create(joinpath(path,name*".nc"),[var,tvar],mode=NC_NETCDF4)
-    NetCDF.putatt(nc,name,Dict("units"=>unit,"long_name"=>long_name))
+    # add missing_value although irrelevant for ncview compatibility
+    NetCDF.putatt(nc,name,Dict("units"=>unit,"long_name"=>long_name,"missing_value"=>-999999f0))
     return nc
 end
 
@@ -107,10 +108,12 @@ function output_nc!(i::Int,
             NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.q != nothing
-            @views q = Float32.(Diag.Vorticity)
+            @views q = Float32.(Diag.Vorticity.q[haloη+1:end-haloη,haloη+1:end-haloη])
             NetCDF.putvar(ncs.q,"q",q,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.ζ != nothing
+            @unpack dvdx,dudy = Diag.Vorticity
+            @unpack f_q = S.grid
             @views ζ = Float32.((dvdx[2:end-1,2:end-1]-dudy[2+ep:end-1,2:end-1])./abs.(f_q))
             NetCDF.putvar(ncs.ζ,"relvort",ζ,start=[1,1,iout],count=[-1,-1,1])
         end