From cdea13d1a28cf7765fbad69d6f6b9f5dc8a73af2 Mon Sep 17 00:00:00 2001 From: Sarah Williamson Date: Tue, 7 May 2024 16:31:11 -0500 Subject: [PATCH 1/5] Adding minor changes that are needed for Enzyme to run, working on understanding why they're needed but for now I just know they are --- src/constants.jl | 2 +- src/default_parameters.jl | 8 ++++---- src/grid.jl | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/constants.jl b/src/constants.jl index 217a44c..73aaf8d 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -26,7 +26,7 @@ function SSPRK3coeff{T}(P::Parameter,Δt_Δ::T) where T return SSPRK3coeff{T}(n,s,kn,mn,Δt_Δn,kna,knb,Δt_Δnc) end -struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat} +mutable struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat} # RUNGE-KUTTA COEFFICIENTS 2nd/3rd/4th order including timestep Δt RKaΔt::Array{Tprog,1} diff --git a/src/default_parameters.jl b/src/default_parameters.jl index 3f54fe3..7df6de2 100644 --- a/src/default_parameters.jl +++ b/src/default_parameters.jl @@ -1,4 +1,4 @@ -@with_kw struct Parameter +@with_kw mutable struct Parameter T=Float32 # number format @@ -8,8 +8,8 @@ # DOMAIN RESOLUTION AND RATIO nx::Int=100 # number of grid cells in x-direction - Lx::Real=4000e3 # length of the domain in x-direction [m] - L_ratio::Real=2 # Domain aspect ratio of Lx/Ly + Lx::Int=4000e3 # length of the domain in x-direction [m] + L_ratio::Int=2 # Domain aspect ratio of Lx/Ly # PHYSICAL CONSTANTS g::Real=0.1 # gravitational acceleration [m/s] @@ -20,7 +20,7 @@ R::Real=6.371e6 # Earth's radius [m] # SCALE - scale::Real=2^6 # multiplicative scale for the momentum equations u,v + scale::Int=2^6 # multiplicative scale for the momentum equations u,v scale_sst::Real=2^15 # multiplicative scale for sst # WIND FORCING OPTIONS diff --git a/src/grid.jl b/src/grid.jl index bab0502..23ae08d 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -2,7 +2,7 @@ # Parameters taken from Parameter struct nx::Int # number of grid cells in x-direction - Lx::Real # length of the domain in x-direction [m] + Lx::Int # length of the domain in x-direction [m] L_ratio::Real # Domain aspect ratio of Lx/Ly bc::String # boundary condition, "periodic" or "nonperiodic" g::Real # gravitational acceleration [m/s] @@ -16,10 +16,10 @@ ω::Real # Earth's angular frequency [s^-1] ϕ::Real # central latitue of the domain (for coriolis) [°] R::Real # Earth's radius [m] - scale::Real # multiplicative scale for momentum equations [1] + scale::Int # multiplicative scale for momentum equations [1] # DOMAIN SIZES - Δ::Real=Lx / nx # grid spacing + Δ::T=Lx / nx # grid spacing ny::Int=Int(round(Lx / L_ratio / Δ)) # number of grid cells in y-direction Ly::Real=ny * Δ # length of domain in y-direction From de2355e55aebeac22587a2e7ed1298632bd255d6 Mon Sep 17 00:00:00 2001 From: Sarah Williamson Date: Tue, 7 May 2024 16:47:10 -0500 Subject: [PATCH 2/5] Fixing the error --- src/grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid.jl b/src/grid.jl index 23ae08d..2fe70db 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -19,7 +19,7 @@ scale::Int # multiplicative scale for momentum equations [1] # DOMAIN SIZES - Δ::T=Lx / nx # grid spacing + Δ::Int=Lx / nx # grid spacing ny::Int=Int(round(Lx / L_ratio / Δ)) # number of grid cells in y-direction Ly::Real=ny * Δ # length of domain in y-direction From cfd1ba2da9c3e0b970524adb42df93440aba3df6 Mon Sep 17 00:00:00 2001 From: Sarah Williamson Date: Wed, 8 May 2024 10:10:20 -0500 Subject: [PATCH 3/5] Making some Ints Float64 --- src/default_parameters.jl | 6 +++--- src/grid.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/default_parameters.jl b/src/default_parameters.jl index 7df6de2..23eea52 100644 --- a/src/default_parameters.jl +++ b/src/default_parameters.jl @@ -8,8 +8,8 @@ # DOMAIN RESOLUTION AND RATIO nx::Int=100 # number of grid cells in x-direction - Lx::Int=4000e3 # length of the domain in x-direction [m] - L_ratio::Int=2 # Domain aspect ratio of Lx/Ly + Lx::Float64=4000e3 # length of the domain in x-direction [m] + L_ratio::Float64=2 # Domain aspect ratio of Lx/Ly # PHYSICAL CONSTANTS g::Real=0.1 # gravitational acceleration [m/s] @@ -20,7 +20,7 @@ R::Real=6.371e6 # Earth's radius [m] # SCALE - scale::Int=2^6 # multiplicative scale for the momentum equations u,v + scale::Float64=2^6 # multiplicative scale for the momentum equations u,v scale_sst::Real=2^15 # multiplicative scale for sst # WIND FORCING OPTIONS diff --git a/src/grid.jl b/src/grid.jl index 2fe70db..3ebbdcd 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -16,7 +16,7 @@ ω::Real # Earth's angular frequency [s^-1] ϕ::Real # central latitue of the domain (for coriolis) [°] R::Real # Earth's radius [m] - scale::Int # multiplicative scale for momentum equations [1] + scale::Float64 # multiplicative scale for momentum equations [1] # DOMAIN SIZES Δ::Int=Lx / nx # grid spacing From 076074c713c005e9ee3456c8ea1895c4e8b63572 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Wed, 8 May 2024 12:32:30 -0400 Subject: [PATCH 4/5] Remove abstract fields from Parameter and Grid --- src/default_parameters.jl | 178 ++++++++++++++++++++------------------ src/grid.jl | 58 ++++++------- 2 files changed, 123 insertions(+), 113 deletions(-) diff --git a/src/default_parameters.jl b/src/default_parameters.jl index 23eea52..774e76b 100644 --- a/src/default_parameters.jl +++ b/src/default_parameters.jl @@ -12,45 +12,45 @@ L_ratio::Float64=2 # Domain aspect ratio of Lx/Ly # PHYSICAL CONSTANTS - g::Real=0.1 # gravitational acceleration [m/s] - H::Real=500. # layer thickness at rest [m] - ρ::Real=1e3 # water density [kg/m^3] - ϕ::Real=45. # central latitude of the domain (for coriolis) [°] - ω::Real=2π/(24*3600) # Earth's angular frequency [s^-1] - R::Real=6.371e6 # Earth's radius [m] + g::Float64 = 0.1 # gravitational acceleration [m/s] + H::Float64 = 500. # layer thickness at rest [m] + ρ::Float64 = 1e3 # water density [kg/m^3] + ϕ::Float64 = 45. # central latitude of the domain (for coriolis) [°] + ω::Float64 = 2π/(24*3600) # Earth's angular frequency [s^-1] + R::Float64 = 6.371e6 # Earth's radius [m] # SCALE scale::Float64=2^6 # multiplicative scale for the momentum equations u,v - scale_sst::Real=2^15 # multiplicative scale for sst + scale_sst::Float64=2^15 # multiplicative scale for sst # WIND FORCING OPTIONS wind_forcing_x::String="shear" # "channel", "double_gyre", "shear","constant" or "none" wind_forcing_y::String="constant" # "channel", "double_gyre", "shear","constant" or "none" - Fx0::Real=0.12 # wind stress strength [Pa] in x-direction - Fy0::Real=0.0 # wind stress strength [Pa] in y-direction + Fx0::Float64=0.12 # wind stress strength [Pa] in x-direction + Fy0::Float64=0.0 # wind stress strength [Pa] in y-direction seasonal_wind_x::Bool=true # Change the wind stress with a sine of frequency ωFx,ωFy seasonal_wind_y::Bool=false # same for y-component - ωFx::Real=2 # frequency [1/year] for x component - ωFy::Real=2 # frequency [1/year] for y component + ωFx::Float64=2 # frequency [1/year] for x component + ωFy::Float64=2 # frequency [1/year] for y component # BOTTOM TOPOGRAPHY OPTIONS topography::String="ridges" # "ridge", "seamount", "flat", "ridges", "bathtub" - topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9] - topo_height::Real=100. # height of seamount [m] - topo_width::Real=300e3 # horizontal scale [m] of the seamount + topo_ridges_positions::Vector{Float64} = [0.05,0.25,0.45,0.9] + topo_height::Float64 = 100. # height of seamount [m] + topo_width::Float64 = 300e3 # horizontal scale [m] of the seamount # SURFACE RELAXATION 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 + t_relax::Float64 = 100. # time scale of the relaxation [days] + η_refh::Float64 = 5. # height difference [m] of the interface relaxation profile + η_refw::Float64 = 50e3 # width [m] of the tangent used for the interface relaxation # SURFACE FORCING surface_forcing::Bool=false # yes? - ωFη::Real=1.0 # frequency [1/year] for surfance forcing - A::Real=3e-5 # Amplitude [m/s] - ϕk::Real=ϕ # Central latitude of Kelvin wave pumping - wk::Real=10e3 # width [m] in y of Gaussian used for surface forcing + ωFη::Float64 = 1.0 # frequency [1/year] for surfance forcing + A::Float64 = 3e-5 # Amplitude [m/s] + ϕk::Float64 = ϕ # Central latitude of Kelvin wave pumping + wk::Float64 = 10e3 # width [m] in y of Gaussian used for surface forcing # TIME STEPPING OPTIONS time_scheme::String="RK" # Runge-Kutta ("RK") or strong-stability preserving RK @@ -58,8 +58,8 @@ RKo::Int=4 # Order of the RK time stepping scheme (2, 3 or 4) RKs::Int=3 # Number of stages for SSPRK2 RKn::Int=5 # n^2 = s = Number of stages for SSPRK3 - cfl::Real=0.9 # CFL number (1.0 recommended for RK4, 0.6 for RK3) - Ndays::Real=200.0 # number of days to integrate for + cfl::Float64 = 0.9 # CFL number (1.0 recommended for RK4, 0.6 for RK3) + Ndays::Float64 = 200.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. # 0 means it is included in every RK4 substep @@ -67,7 +67,7 @@ # BOUNDARY CONDITION OPTIONS bc::String="periodic" # "periodic" or anything else for nonperiodic - α::Real=2. # lateral boundary condition parameter + α::Float64 = 2 # lateral boundary condition parameter # 0 free-slip, 0<α<2 partial-slip, 2 no-slip # PARAMETERS FOR ADJOINT METHOD @@ -85,13 +85,13 @@ # BOTTOM FRICTION OPTIONS bottom_drag::String="none" # "linear", "quadratic" or "none" - cD::Real=1e-5 # bottom drag coefficient [dimensionless] for quadratic - τD::Real=300. # bottom drag coefficient [days] for linear + cD::Float64 = 1e-5 # bottom drag coefficient [dimensionless] for quadratic + τD::Float64 = 300. # bottom drag coefficient [days] for linear # DIFFUSION OPTIONS diffusion::String="constant" # "Smagorinsky" or "constant", biharmonic in both cases - νB::Real=500.0 # [m^2/s] scaling constant for constant biharmonic diffusion - cSmag::Real=0.15 # Smagorinsky coefficient [dimensionless] + νB::Float64 = 500.0 # [m^2/s] scaling constant for constant biharmonic diffusion + cSmag::Float64 = 0.15 # Smagorinsky coefficient [dimensionless] # TRACER ADVECTION tracer_advection::Bool=true # yes? @@ -100,21 +100,21 @@ sst_initial::String="waves" # "west", "south", "linear", "waves","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.2 # Velocity scale [m/s] for tracer advection - SSTmax::Real=1. # tracer (sea surface temperature) max for initial conditions - SSTmin::Real=-1. # tracer (sea surface temperature) min for initial conditions - τSST::Real=100 # tracer restoring time scale [days] - jSST::Real=365 # tracer consumption [days] - SSTw::Real=5e5 # width [m] of the tangent used for the IC and interface relaxation - SSTϕ::Real=0.5 # latitude/longitude fraction ∈ [0,1] of sst edge - SSTwaves_ny::Real=4 # wave crests/troughs in y - SSTwaves_nx::Real=SSTwaves_ny*L_ratio # wave crests/troughs in x - SSTwaves_p::Real=1/2 # power for rectangles (p<1)/smootheness(p>=1) of waves + Uadv::Float64 = 0.2 # Velocity scale [m/s] for tracer advection + SSTmax::Float64 = 1. # tracer (sea surface temperature) max for initial conditions + SSTmin::Float64 = -1. # tracer (sea surface temperature) min for initial conditions + τSST::Float64 = 100 # tracer restoring time scale [days] + jSST::Float64 = 365 # tracer consumption [days] + SSTw::Float64 = 5e5 # width [m] of the tangent used for the IC and interface relaxation + SSTϕ::Float64 = 0.5 # latitude/longitude fraction ∈ [0,1] of sst edge + SSTwaves_ny::Float64 = 4 # wave crests/troughs in y + SSTwaves_nx::Float64 = SSTwaves_ny*L_ratio # wave crests/troughs in x + SSTwaves_p::Float64 = 1/2 # power for rectangles (p<1)/smootheness(p>=1) of waves # OUTPUT OPTIONS output::Bool=false # netcdf output? output_vars::Array{String,1}=["u","v","η","sst"] # which variables to output? "du","dv","dη","H","ζ" also allowed - output_dt::Real=24 # output time step [hours] + output_dt::Float64 = 24 # output time step [hours] outpath::String=pwd() # path to output folder compression_level::Int=3 # compression level return_time::Bool=false # return time of simulation of progn vars? @@ -164,7 +164,6 @@ end """ Creates a Parameter struct with following options and default values - T=Float32 # number format Tprog=T # number format for prognostic variables @@ -173,80 +172,90 @@ Creates a Parameter struct with following options and default values # DOMAIN RESOLUTION AND RATIO nx::Int=100 # number of grid cells in x-direction - Lx::Real=4000e3 # length of the domain in x-direction [m] - L_ratio::Real=2 # Domain aspect ratio of Lx/Ly + Lx::Float64=4000e3 # length of the domain in x-direction [m] + L_ratio::Float64=2 # Domain aspect ratio of Lx/Ly # PHYSICAL CONSTANTS - g::Real=0.1 # gravitational acceleration [m/s] - H::Real=500. # layer thickness at rest [m] - ρ::Real=1e3 # water density [kg/m^3] - ϕ::Real=45. # central latitude of the domain (for coriolis) [°] - ω::Real=2π/(24*3600) # Earth's angular frequency [s^-1] - R::Real=6.371e6 # Earth's radius [m] + g::Float64 = 0.1 # gravitational acceleration [m/s] + H::Float64 = 500. # layer thickness at rest [m] + ρ::Float64 = 1e3 # water density [kg/m^3] + ϕ::Float64 = 45. # central latitude of the domain (for coriolis) [°] + ω::Float64 = 2π/(24*3600) # Earth's angular frequency [s^-1] + R::Float64 = 6.371e6 # Earth's radius [m] # SCALE - scale::Real=2^6 # multiplicative scale for the momentum equations u,v - scale_sst::Real=2^15 # multiplicative scale for sst + scale::Float64=2^6 # multiplicative scale for the momentum equations u,v + scale_sst::Float64=2^15 # multiplicative scale for sst # WIND FORCING OPTIONS wind_forcing_x::String="shear" # "channel", "double_gyre", "shear","constant" or "none" wind_forcing_y::String="constant" # "channel", "double_gyre", "shear","constant" or "none" - Fx0::Real=0.12 # wind stress strength [Pa] in x-direction - Fy0::Real=0.0 # wind stress strength [Pa] in y-direction + Fx0::Float64=0.12 # wind stress strength [Pa] in x-direction + Fy0::Float64=0.0 # wind stress strength [Pa] in y-direction seasonal_wind_x::Bool=true # Change the wind stress with a sine of frequency ωFx,ωFy seasonal_wind_y::Bool=false # same for y-component - ωFx::Real=2 # frequency [1/year] for x component - ωFy::Real=2 # frequency [1/year] for y component + ωFx::Float64=2 # frequency [1/year] for x component + ωFy::Float64=2 # frequency [1/year] for y component # BOTTOM TOPOGRAPHY OPTIONS topography::String="ridges" # "ridge", "seamount", "flat", "ridges", "bathtub" - topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9] - topo_height::Real=100. # height of seamount [m] - topo_width::Real=300e3 # horizontal scale [m] of the seamount + topo_ridges_positions::Vector{Float64} = [0.05,0.25,0.45,0.9] + topo_height::Float64 = 100. # height of seamount [m] + topo_width::Float64 = 300e3 # horizontal scale [m] of the seamount # SURFACE RELAXATION 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 + t_relax::Float64 = 100. # time scale of the relaxation [days] + η_refh::Float64 = 5. # height difference [m] of the interface relaxation profile + η_refw::Float64 = 50e3 # width [m] of the tangent used for the interface relaxation # SURFACE FORCING surface_forcing::Bool=false # yes? - ωFη::Real=1.0 # frequency [1/year] for surfance forcing - A::Real=3e-5 # Amplitude [m/s] - ϕk::Real=ϕ # Central latitude of Kelvin wave pumping - wk::Real=10e3 # width [m] in y of Gaussian used for surface forcing + ωFη::Float64 = 1.0 # frequency [1/year] for surfance forcing + A::Float64 = 3e-5 # Amplitude [m/s] + ϕk::Float64 = ϕ # Central latitude of Kelvin wave pumping + wk::Float64 = 10e3 # width [m] in y of Gaussian used for surface forcing # TIME STEPPING OPTIONS - time_scheme::String="SSPRK3" # Runge-Kutta ("RK") or strong-stability preserving RK + time_scheme::String="RK" # Runge-Kutta ("RK") or strong-stability preserving RK # "SSPRK2","SSPRK3","4SSPRK3" RKo::Int=4 # Order of the RK time stepping scheme (2, 3 or 4) RKs::Int=3 # Number of stages for SSPRK2 RKn::Int=5 # n^2 = s = Number of stages for SSPRK3 - cfl::Real=4.0 # CFL number (1.0 recommended for RK4, 0.6 for RK3) - Ndays::Real=200.0 # number of days to integrate for + cfl::Float64 = 0.9 # CFL number (1.0 recommended for RK4, 0.6 for RK3) + Ndays::Float64 = 200.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. # 0 means it is included in every RK4 substep + compensated::Bool=false # Compensated summation in the time integration? # BOUNDARY CONDITION OPTIONS bc::String="periodic" # "periodic" or anything else for nonperiodic - α::Real=2. # lateral boundary condition parameter + α::Float64 = 2 # lateral boundary condition parameter # 0 free-slip, 0<α<2 partial-slip, 2 no-slip + # PARAMETERS FOR ADJOINT METHOD + data_steps::StepRange{Int,Int} = 0:1:0 # Timesteps where data exists + data::Array{Float32, 1} = [0.] # model data + J::Float64 = 0. # Placeholder for cost function evaluation + j::Int = 1 # For keeping track of the entry in data + + # CHECKPOINTING VARIABLES + i::Int = 0 # Placeholder for current timestep, needed for Checkpointing.jl + # MOMENTUM ADVECTION OPTIONS adv_scheme::String="ArakawaHsu" # "Sadourny" or "ArakawaHsu" dynamics::String="nonlinear" # "linear" or "nonlinear" # BOTTOM FRICTION OPTIONS bottom_drag::String="none" # "linear", "quadratic" or "none" - cD::Real=1e-5 # bottom drag coefficient [dimensionless] for quadratic - τD::Real=300. # bottom drag coefficient [days] for linear + cD::Float64 = 1e-5 # bottom drag coefficient [dimensionless] for quadratic + τD::Float64 = 300. # bottom drag coefficient [days] for linear # DIFFUSION OPTIONS diffusion::String="constant" # "Smagorinsky" or "constant", biharmonic in both cases - νB::Real=500.0 # [m^2/s] scaling constant for constant biharmonic diffusion - cSmag::Real=0.15 # Smagorinsky coefficient [dimensionless] + νB::Float64 = 500.0 # [m^2/s] scaling constant for constant biharmonic diffusion + cSmag::Float64 = 0.15 # Smagorinsky coefficient [dimensionless] # TRACER ADVECTION tracer_advection::Bool=true # yes? @@ -255,23 +264,24 @@ Creates a Parameter struct with following options and default values sst_initial::String="waves" # "west", "south", "linear", "waves","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.2 # Velocity scale [m/s] for tracer advection - SSTmax::Real=1. # tracer (sea surface temperature) max for initial conditions - SSTmin::Real=-1. # tracer (sea surface temperature) min for initial conditions - τSST::Real=100 # tracer restoring time scale [days] - jSST::Real=365 # tracer consumption [days] - SSTw::Real=5e5 # width [m] of the tangent used for the IC and interface relaxation - SSTϕ::Real=0.5 # latitude/longitude fraction ∈ [0,1] of sst edge - SSTwaves_ny::Real=4 # wave crests/troughs in y - SSTwaves_nx::Real=SSTwaves_ny*L_ratio # wave crests/troughs in x - SSTwaves_p::Real=1/2 # power for rectangles (p<1)/smootheness(p>=1) of waves + Uadv::Float64 = 0.2 # Velocity scale [m/s] for tracer advection + SSTmax::Float64 = 1. # tracer (sea surface temperature) max for initial conditions + SSTmin::Float64 = -1. # tracer (sea surface temperature) min for initial conditions + τSST::Float64 = 100 # tracer restoring time scale [days] + jSST::Float64 = 365 # tracer consumption [days] + SSTw::Float64 = 5e5 # width [m] of the tangent used for the IC and interface relaxation + SSTϕ::Float64 = 0.5 # latitude/longitude fraction ∈ [0,1] of sst edge + SSTwaves_ny::Float64 = 4 # wave crests/troughs in y + SSTwaves_nx::Float64 = SSTwaves_ny*L_ratio # wave crests/troughs in x + SSTwaves_p::Float64 = 1/2 # power for rectangles (p<1)/smootheness(p>=1) of waves # OUTPUT OPTIONS output::Bool=false # netcdf output? - output_vars::Array{String,1}=["u","v","η","sst"] # which variables to output? "du","dv","dη" also allowed - output_dt::Real=24 # output time step [hours] + output_vars::Array{String,1}=["u","v","η","sst"] # which variables to output? "du","dv","dη","H","ζ" also allowed + output_dt::Float64 = 24 # output time step [hours] outpath::String=pwd() # path to output folder compression_level::Int=3 # compression level + return_time::Bool=false # return time of simulation of progn vars? # INITIAL CONDITIONS initial_cond::String="rest" # "rest" or "ncfile" for restart from file diff --git a/src/grid.jl b/src/grid.jl index 3ebbdcd..162f8f1 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -3,25 +3,25 @@ # Parameters taken from Parameter struct nx::Int # number of grid cells in x-direction Lx::Int # length of the domain in x-direction [m] - L_ratio::Real # Domain aspect ratio of Lx/Ly + L_ratio::Float64 # Domain aspect ratio of Lx/Ly bc::String # boundary condition, "periodic" or "nonperiodic" - g::Real # gravitational acceleration [m/s] - H::Real # layer thickness at rest [m] - cfl::Real # CFL number (1.0 recommended for RK4, 0.6 for RK3) - Ndays::Real # number of days to integrate for + g::Float64 # gravitational acceleration [m/s] + H::Float64 # layer thickness at rest [m] + cfl::Float64 # CFL number (1.0 recommended for RK4, 0.6 for RK3) + Ndays::Float64 # number of days to integrate for nstep_diff::Int # diffusive terms every nstep_diff time steps. nstep_advcor::Int # nonlinear terms every nstep_advcor time steps. - Uadv::Real # Velocity scale [m/s] for tracer advection - output_dt::Real # output time step [hours] - ω::Real # Earth's angular frequency [s^-1] - ϕ::Real # central latitue of the domain (for coriolis) [°] - R::Real # Earth's radius [m] + Uadv::Float64 # Velocity scale [m/s] for tracer advection + output_dt::Float64 # output time step [hours] + ω::Float64 # Earth's angular frequency [s^-1] + ϕ::Float64 # central latitue of the domain (for coriolis) [°] + R::Float64 # Earth's radius [m] scale::Float64 # multiplicative scale for momentum equations [1] # DOMAIN SIZES - Δ::Int=Lx / nx # grid spacing + Δ::Int=Lx / nx # grid spacing ny::Int=Int(round(Lx / L_ratio / Δ)) # number of grid cells in y-direction - Ly::Real=ny * Δ # length of domain in y-direction + Ly::Float64=ny * Δ # length of domain in y-direction # NUMBER OF GRID POINTS nux::Int = if (bc == "periodic") nx else nx-1 end # u-grid in x-direction @@ -38,14 +38,14 @@ nq::Int = nqx*nqy # q-grid # GRID VECTORS - x_T::AbstractVector = Δ*Array(1:nx) .- Δ/2 - y_T::AbstractVector = Δ*Array(1:ny) .- Δ/2 - x_u::AbstractVector = if (bc == "periodic") Δ*Array(0:nx-1) else Δ*Array(1:nx-1) end - y_u::AbstractVector = y_T - x_v::AbstractVector = x_T - y_v::AbstractVector = Δ*Array(1:ny-1) - x_q::AbstractVector = if bc == "periodic" x_u else Δ*Array(1:nx+1) .- Δ end - y_q::AbstractVector = Δ*Array(1:ny+1) .- Δ + x_T::Vector{Float64} = Δ*Array(1:nx) .- Δ/2 + y_T::Vector{Float64} = Δ*Array(1:ny) .- Δ/2 + x_u::Vector{Float64} = if (bc == "periodic") Δ*Array(0:nx-1) else Δ*Array(1:nx-1) end + y_u::Vector{Float64} = y_T + x_v::Vector{Float64} = x_T + y_v::Vector{Float64} = Δ*Array(1:ny-1) + x_q::Vector{Float64} = if bc == "periodic" x_u else Δ*Array(1:nx+1) .- Δ end + y_q::Vector{Float64} = Δ*Array(1:ny+1) .- Δ # HALO SIZES halo::Int=2 # halo size for u,v (Biharmonic stencil requires 2) @@ -57,17 +57,17 @@ ep::Int = if bc == "periodic" 1 else 0 end # is there a u-point on the left edge? # GRID VECTORS WITH HALO - x_T_halo::AbstractVector = Δ*Array(0:nx+1) .- Δ/2 - y_T_halo::AbstractVector = Δ*Array(0:ny+1) .- Δ/2 - x_u_halo::AbstractVector = if (bc == "periodic") Δ*Array(-2:nx+1) else Δ*Array(-1:nx+1) end - y_u_halo::AbstractVector = Δ*Array(-1:ny+2) .- Δ/2 - x_v_halo::AbstractVector = Δ*Array(-1:nx+2) .- Δ/2 - y_v_halo::AbstractVector = Δ*Array(-1:ny+1) - x_q_halo::AbstractVector = if bc == "periodic" x_u_halo else Δ*Array(-1:nx+3) .- Δ end - y_q_halo::AbstractVector = Δ*Array(-1:ny+3) .- Δ + x_T_halo::Vector{Float64} = Δ*Array(0:nx+1) .- Δ/2 + y_T_halo::Vector{Float64} = Δ*Array(0:ny+1) .- Δ/2 + x_u_halo::Vector{Float64} = if (bc == "periodic") Δ*Array(-2:nx+1) else Δ*Array(-1:nx+1) end + y_u_halo::Vector{Float64} = Δ*Array(-1:ny+2) .- Δ/2 + x_v_halo::Vector{Float64} = Δ*Array(-1:nx+2) .- Δ/2 + y_v_halo::Vector{Float64} = Δ*Array(-1:ny+1) + x_q_halo::Vector{Float64} = if bc == "periodic" x_u_halo else Δ*Array(-1:nx+3) .- Δ end + y_q_halo::Vector{Float64} = Δ*Array(-1:ny+3) .- Δ # TIME STEPS - c::Real = √(g*H) # shallow water gravity wave speed + c::Float64 = √(g*H) # shallow water gravity wave speed dtint::Int = Int(floor(cfl*Δ/c)) # dt converted to Int nt::Int = Int(ceil(Ndays*3600*24/dtint)) # number of time steps to integrate dt::T = T(dtint) # time step [s] From ffdc31e874833d8294d9180082c058bffc2e241f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Wed, 8 May 2024 13:45:34 -0400 Subject: [PATCH 5/5] bump version to 0.5.2 --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index bf5e042..e79799a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ShallowWaters" uuid = "56019723-2d87-4a65-81ff-59d5d8913e3c" -authors = ["Milan Kloewer"] -version = "0.5.1" +authors = ["Milan Kloewer and ShallowWaters.jl contributors"] +version = "0.5.2" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"