Merge pull request #109 from milankl/output

Output

src/DefaultParameters.jl

@@ -23,7 +23,7 @@
 
     # BOTTOM TOPOGRAPHY OPTIONS
     topography::String="ridge"          # "ridge", "seamount", "flat", "ridges", "bathtub"
-    topo_height::Real=10.               # height of seamount [m]
+    topo_height::Real=50.               # height of seamount [m]
     topo_width::Real=300e3              # horizontal scale [m] of the seamount
 
     # SURFACE RELAXATION
@@ -61,54 +61,66 @@
 
     # DIFFUSION OPTIONS
     diffusion::String="Smagorinsky"     # "Smagorinsky" or "Constant", biharmonic in both cases
-    νB::Real=500.0                      # [m^2/s] scaling constant for Constant biharmonic diffusion
+    νB::Real=500.0                      # [m^2/s] scaling constant for constant biharmonic diffusion
     cSmag::Real=0.15                    # Smagorinsky coefficient [dimensionless]
 
     # TRACER ADVECTION
-    tracer_advection::Bool=false        # yes?
+    tracer_advection::Bool=true         # yes?
     tracer_relaxation::Bool=false       # yes?
     tracer_consumption::Bool=false      # yes?
     tracer_pumping::Bool=false          # yes?
-    injection_region::String="west"     # "west", "south", "rect" or flat
+    injection_region::String="west"     # "west" or "south"
     sst_initial::String="south"         # same here
     sstrestart::Bool=true               # start from previous sst file
-    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]
     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
-    SSTw::Real=1000e3                   # width [m] of the tangent used for the IC and interface relaxation
+    SSTw::Real=50e3                     # width [m] of the tangent used for the IC and interface relaxation
     SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
 
     # OUTPUT OPTIONS
-    output::Bool=false                  # for nc output
-    output_tend::Bool=false             # ouput for tendencies as well?
-    output_diagn::Bool=false            # output diagnostic variables as well?
-
-                                        # which prognostic variables to output?
-    output_progn_vars::Array{String,1}=["u","v","eta","sst"]
-                                        # which tendencies to output?
-    output_tend_vars::Array{String,1}=["du","dv","deta","Bu","Bv","LLu1","LLu2","LLv1","LLv2",
-                                "qhv","qhu","dpdx","dpdy","dUdx","dVdy"]
-                                        # which diagnostic variables to output?
-    output_diagn_vars::Array{String,1}=["q","p","dudx","dvdy","dudy","dvdx","Lu",
-                                "Lv","xd","yd"]
-
+    output::Bool=false                  # netcdf output?
+    output_vars::Array{String,1}=["u","v","eta","sst","q","ζ"]  # which variables to output?
     output_dt::Real=6                   # output time step [hours]
-    outpath::String="data/"             # path to output folder
+    outpath::String=pwd()               # path to output folder
 
     # INITIAL CONDITIONS
     initial_cond::String="rest"         # "rest" or "ncfile" for restart from file
-    initpath::String="data/"            # folder where to pick the restart files from
+    initpath::String=outpath            # folder where to pick the restart files from
     init_run_id::Int=0                  # run id for restart from run number
 
     # ASSERT - CHECK THAT THE INPUT PARAMETERS MAKE SENSE
-    @assert all((nx,Lx,L_ratio) .> 0.)
-    @assert all((g,H,ρ,ω,R) .> 0.)
-    @assert ϕ <= 90.0 && ϕ >= -90.0
-    #TODO more of that
-
+    @assert all((nx,Lx,L_ratio) .> 0.)  "nx, Lx, L_ratio have to be >0"
+    @assert all((g,H,ρ,ω,R) .> 0.)      "g,H,ρ,ω,R have to be >0"
+    @assert ϕ <= 90.0 && ϕ >= -90.0     "ϕ has to be in (-90,90), $ϕ given."
+    @assert wind_forcing_x in ["channel","double_gyre","shear","constant","none"] "Wind forcing '$wind_forcing_x' unsupported"
+    @assert wind_forcing_y in ["channel","double_gyre","shear","constant","none"] "Wind forcing '$wind_forcing_y' unsupported"
+    @assert topography in ["ridge","seamount","flat","ridges"] "Topography '$topography' unsupported"
+    @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 ωyr > 0.0           "ωyr has to be >0, $ωyr given."
+    @assert RKo in [3,4]        "RKo has to be 3 or 4, $RKo given."
+    @assert Ndays > 0.0         "Ndays has to be >0, $Ndyas 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."
+    @assert bc in ["periodic","nonperiodic"]    "boundary condition '$bc' unsupported."
+    @assert α >= 0.0 && α <= 2.0    "Tangential boundary condition α has to be in (0,2), $α given."
+    @assert adv_scheme in ["Sadourny","ArakawaHsu"] "Advection scheme '$adv_scheme' unsupported"
+    @assert dynamics in ["linear","nonlinear"]  "Dynamics '$dynamics' unsupoorted."
+    @assert bottom_drag in ["quadratic","linear","none"] "Bottom drag '$bottom_drag' unsupported."
+    @assert cD >= 0.0    "Bottom drag coefficient cD has to be >=0, $cD given."
+    @assert τD >= 0.0    "Bottom drag coefficient τD has to be >=0, $τD given."
+    @assert diffusion in ["Smagorinsky", "constant"] "Diffusion '$diffusion' unsupoorted."
+    @assert νB > 0.0     "Diffusion scaling constant νB has to be > 0, $νB given."
+    @assert cSmag > 0.0  "Smagorinsky coefficient cSmag has to be >0, $cSmag given."
+    @assert injection_region in ["west","south"] "Injection region '$injection_region' unsupported."
+    @assert Uadv > 0.0   "Advection velocity scale Uadv has to be >0, $Uadv given."
+    @assert output_dt > 0   "Output time step has to be >0, $output_dt given."
+    @assert initial_cond in ["rest", "ncfile"] "Initial conditions '$initial_cond' unsupported."
 end

src/Diffusion.jl

@@ -8,6 +8,8 @@ function diffusion!(    u::AbstractMatrix,
         diffusion_constant!(u,v,Diag,S)
     elseif S.parameters.diffusion == "Smagorinsky"
         diffusion_smagorinsky!(u,v,Diag,S)
+    else
+        throw(error("Diffusion scheme $(S.parameters.diffusion) is unsupported."))
     end
 end
 
@@ -168,55 +170,6 @@ function viscous_tensor_smagorinsky!(Diag::DiagnosticVars)
     end
 end
 
-"""Biharmonic stress tensor times Smagorinsky coefficient
-νSmag * ∇∇² ⃗u = (S11, S12; S21, S22)."""
-function viscous_tensor_smagorinsky!(Diag::DiagnosticVars)
-
-    @unpack dLudx,dLudy,dLvdx,dLvdy = Diag.Laplace
-    @unpack νSmag,νSmag_q,S11,S12,S21,S22 = Diag.Smagorinsky
-    @unpack ep = Diag.Laplace
-
-    m,n = size(S11)
-    @boundscheck (m+2-ep,n) == size(νSmag) || throw(BoundsError())
-    @boundscheck (m,n) == size(dLudx) || throw(BoundsError())
-
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            S11[i,j] = νSmag[i+1-ep,j] * dLudx[i,j]
-        end
-    end
-
-    m,n = size(S12)
-    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
-    @boundscheck (m+ep,n) == size(dLudy) || throw(BoundsError())
-
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            S12[i,j] = νSmag_q[i,j] * dLudy[i+ep,j]
-        end
-    end
-
-    m,n = size(S21)
-    @boundscheck (m,n) == size(νSmag_q) || throw(BoundsError())
-    @boundscheck (m,n) == size(dLvdx) || throw(BoundsError())
-
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            S21[i,j] = νSmag_q[i,j] * dLvdx[i,j]
-        end
-    end
-
-    m,n = size(S22)
-    @boundscheck (m,n+2) == size(νSmag) || throw(BoundsError())
-    @boundscheck (m,n) == size(dLvdy) || throw(BoundsError())
-
-    @inbounds for j ∈ 1:n
-        for i ∈ 1:m
-            S22[i,j] = νSmag[i,j+1] * dLvdy[i,j]
-        end
-    end
-end
-
 """Biharmonic stress tensor times constant viscosity coefficient
 νB * ∇∇² ⃗u = (S11, S12; S21, S22)"""
 function viscous_tensor_constant!(  Diag::DiagnosticVars,

src/Feedback.jl

@@ -6,6 +6,8 @@
     i::Int=0                                # time step increment
     nt::Int                                 # number of time steps
     nout::Int                               # number of time steps with output
+    run_id::Int                             # run identification number
+    runpath::String                         # output path plus run????/
 end
 
 
@@ -32,8 +34,12 @@ function readable_secs(secs::Real)
 end
 
 """Estimates the total time the model integration will take."""
-function duration_estimate(i,t,nt,progrtxt,P::Parameter)
-    time_per_step = (time()-t) / (i-nadvstep)
+function duration_estimate(feedback::Feedback,S::ModelSetup)
+
+    @unpack t0,i,nt,output,progress_txt = feedback
+    @unpack nadvstep = S.grid
+
+    time_per_step = (time()-t0) / (i-nadvstep)
     time_total = Int(round(time_per_step*nt))
     time_to_go = Int(round(time_per_step*(nt-i)))
 
@@ -43,9 +49,9 @@ function duration_estimate(i,t,nt,progrtxt,P::Parameter)
     println(s1)     # print inline
     println(s2)
     if output == 1  # print in txt
-        write(progrtxt,"\n"*s1*"\n")
-        write(progrtxt,s2*"\n")
-        flush(progrtxt)
+        write(progress_txt,"\n"*s1*"\n")
+        write(progress_txt,s2*"\n")
+        flush(progress_txt)
     end
 end
 
@@ -69,7 +75,13 @@ function feedback_init(S::ModelSetup)
     @unpack nt,nout = S.grid
 
     if output
-        txt = open(runpath*"progress.txt","w")
+
+        @unpack Ndays,initial_cond,init_run_id,T,bc,α = S.parameters
+        @unpack nx,ny,Δ = S.grid
+
+        run_id,runpath = get_run_id_path(S)
+
+        txt = open(joinpath(runpath,"progress.txt"),"w")
         s = "Starting Juls run $run_id on "*Dates.format(now(),Dates.RFC1123Format)
         println(s)
         write(txt,s*"\n")
@@ -80,29 +92,31 @@ function feedback_init(S::ModelSetup)
         else
             write(txt,"last time step of run $init_run_id.\n")
         end
-        write(txt,"Boundary conditions are $bc_x with lbc=$lbc.\n")
-        write(txt,"Numtype is "*string(Numtype)*".\n")
-        write(txt,"Time steps are (Lin,Diff,Advcor,Lagr,Output)\n")
-        write(txt,"$dtint, $(dtint*nstep_diff), $(dtint*nstep_advcor), $dtadvint, $(output_dt*3600)\n")
+        write(txt,"Boundary conditions are $bc with lbc=$α.\n")
+        write(txt,"Number format is "*string(T)*".\n")
         write(txt,"\nAll data will be stored in $runpath\n")
     else
         println("Starting Juls on "*Dates.format(now(),Dates.RFC1123Format)*" without output.")
         txt = nothing
+        run_id = -1
+        runpath = ""
     end
 
-    return Feedback(progress_txt=txt,output=output,nt=nt,nout=nout)
+    return Feedback(progress_txt=txt,output=output,nt=nt,nout=nout,run_id=run_id,runpath=runpath)
 end
 
 """Feedback function that calls duration estimate, nan_detection and progress."""
-function feedback!(Prog::PrognosticVars,feedback::Feedback)
+function feedback!(Prog::PrognosticVars,feedback::Feedback,S::ModelSetup)
+
+    @unpack nadvstep = S.grid
+    @unpack i = feedback
 
-    # if i == nadvstep # measure time after tracer advection executed once
-    #     t = time()
-    # elseif i == min(2*nadvstep,nadvstep+50)
-    #     # after the tracer advection executed twice or at least 50 steps
-    #     duration_estimate(i,t,nt,progrtxt)
-    # end
-    #
+    if i == nadvstep # measure time after tracer advection executed once
+        feedback.t0 = time()
+    elseif i == min(2*nadvstep,nadvstep+50)
+        # after the tracer advection executed twice or at least 50 steps
+        duration_estimate(feedback,S)
+    end
 
     @unpack i, nans_detected, nout, output = feedback
 
@@ -131,8 +145,8 @@ function feedback_end!(feedback::Feedback)
     s = " Integration done in "*readable_secs(time()-t0)*"."
     println(s)
     if output
-        write(progrtxt,"\n"*s[2:end]*"\n")  # close txt file with last output
-        flush(progrtxt)
+        write(progress_txt,"\n"*s[2:end]*"\n")  # close txt file with last output
+        flush(progress_txt)
     end
 end
 

src/Grid.jl

@@ -79,6 +79,8 @@
     dtadvint::Int = nadvstep*dtint                          # advection time step [s]
     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 
 
     # N TIME STEPS FOR OUTPUT
     nout::Int = Int(floor(output_dt*3600/dtint))            # output every n time steps

src/Juls.jl

@@ -2,9 +2,7 @@ module Juls
 
 export RunJuls, Parameter
 
-#using NetCDF
-#using Dates, FileIO, Printf, Parameters
-using Parameters, Printf, Dates
+using NetCDF, Parameters, Printf, Dates
 
 include("DefaultParameters.jl")
 include("Grid.jl")
@@ -23,10 +21,10 @@ include("PVadvection.jl")
 include("Continuity.jl")
 include("Bottomdrag.jl")
 include("Diffusion.jl")
-#include("TracerAdvection.jl")
+include("TracerAdvection.jl")
 
 include("Feedback.jl")
-#include("Output.jl")
+include("Output.jl")
 include("RunJuls.jl")
 
 end