Rename DNS - > DNSModel, TLDNS_simulation_setup + update README.md + change default behaviour of overwrite_existing

Project.toml

@@ -1,7 +1,7 @@
 name = "TwoLayerDirectNumericalShenanigans"
 uuid = "40aaee9f-3595-48be-b36c-f1067009652f"
 authors = ["Josef Bisits <[email protected]>"]
-version = "0.4.5"
+version = "0.5.0"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

README.md

@@ -2,7 +2,7 @@
 
 [![Build Status](https://github.com/jbisits/TwoLayerDirectNumericalShenanigans.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/jbisits/TwoLayerDirectNumericalShenanigans.jl/actions/workflows/CI.yml?query=branch%3Amain)
 
-A Julia package (unregistered) to explore the cabbeling instability using Direct Numerical Simulation experiments built with [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl).
+A Julia package (unregistered) to explore two layer systems using Direct Numerical Simulation experiments built with [Oceananigans.jl](https://github.com/CliMA/Oceananigans.jl).
 
 To add the package (assuming julia is already installed):
 

examples/twolayer_example.jl

@@ -6,32 +6,31 @@ diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
 resolution = (Nx = 10, Ny = 10, Nz = 100)
 
 ## Setup the model
-model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
-            reference_density = REFERENCE_DENSITY)
+model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+                 reference_density = REFERENCE_DENSITY)
 
 ## set initial conditions
 T₀ᵘ = -1.5
 S₀ᵘ = (stable = 34.551, cabbeling = 34.58, unstable = 34.59)
 cabbeling = CabbelingUpperLayerInitialConditions(S₀ᵘ.cabbeling, T₀ᵘ)
 initial_conditions = TwoLayerInitialConditions(cabbeling)
 transition_depth = find_depth(model, INTERFACE_LOCATION)
-profile_function = StepChange(transition_depth)#HyperbolicTangent(INTERFACE_LOCATION, 1.0)
+profile_function = StepChange(transition_depth)
 tracer_perturbation_depth = find_depth(model, INTERFACE_LOCATION / 2)
 tracer_perturbation = SalinityGaussianProfile(tracer_perturbation_depth, 0.0, 1.5)
 noise_depth = find_depth(model, INTERFACE_LOCATION)
 initial_noise = SalinityNoise(noise_depth, 1e-2)
 
-dns = TwoLayerDNS(model, profile_function, initial_conditions; initial_noise)
+tldns = TwoLayerDNS(model, profile_function, initial_conditions; initial_noise)
 
-set_two_layer_initial_conditions!(dns)
+set_two_layer_initial_conditions!(tldns)
 
 ## build the simulation
 Δt = 1e-4
 stop_time = 60
 save_schedule = 0.5 # seconds
-output_path = joinpath(@__DIR__, "outputs/")
-simulation = DNS_simulation_setup(dns, Δt, stop_time, save_schedule,
-                                  overwrite_existing = true; output_path)
+output_path = joinpath(@__DIR__, "outputs")
+simulation = TLDNS_simulation_setup(tldns, Δt, stop_time, save_schedule; output_path)
 
 ## Run the simulation
 run!(simulation)

src/TwoLayerDirectNumericalShenanigans.jl

@@ -28,7 +28,7 @@ abstract type AbstractProfileFunction end
 "Abstract supertype for dns initial conditions."
 abstract type AbstractInitialConditions end
 
-export TwoLayerDNS, DNS, DNS_simulation_setup
+export TwoLayerDNS, DNSModel, TLDNS_simulation_setup
 
 export
     StableUpperLayerInitialConditions,

src/twolayerdns.jl

@@ -33,9 +33,10 @@ Base.iterate(tldns::TwoLayerDNS, state = 1) =
     state > length(fieldnames(TwoLayerDNS)) ? nothing :
                                             (getfield(tldns, state), state + 1)
 """
-    function DNS(architecture, domain_extent::NamedTuple, resolution::NamedTuple,
-                 diffusivities::NamedTuple)
-Setup a Direct Numerical Simulation on `architecture` (`CPU()` or `GPU()`) over the
+    function DNSModel(architecture, domain_extent::NamedTuple, resolution::NamedTuple,
+                      diffusivities::NamedTuple)
+Setup a Direct Numerical Simulation [`Model`](https://clima.github.io/OceananigansDocumentation/dev/appendix/library/#Oceananigans.Models.NonhydrostaticModels.NonhydrostaticModel-Tuple{})
+on `architecture` (`CPU()` or `GPU()`) over the
 `domain_extent` with `resolution` and scalar diffusivities for momentum
 (kinematic viscosity `ν`) and the temperature and salinity tracers (`κₜ` and `κₛ`).
 To evolve the temperature and salinity tracers, the [polynomial approximation to the
@@ -66,15 +67,15 @@ the rate `stretching`, if `false` uniform grid spacing is used;
 - `refinement = 1.2` spacing near the surface in the `z` dimension;
 - `stretching = 100` rate of stretching at the bottom of grid in the `z` dimension.
 """
-function DNS(architecture, domain_extent::NamedTuple, resolution::NamedTuple,
-             diffusivities::NamedTuple;
-             linear_eos = false,
-             α = 1.67e-4,
-             β = 7.80e-4,
-             reference_density = nothing,
-             zgrid_stretching = false,
-             refinement = 1.05,
-             stretching = 40)
+function DNSModel(architecture, domain_extent::NamedTuple, resolution::NamedTuple,
+                  diffusivities::NamedTuple;
+                  linear_eos = false,
+                  α = 1.67e-4,
+                  β = 7.80e-4,
+                  reference_density = nothing,
+                  zgrid_stretching = false,
+                  refinement = 1.05,
+                  stretching = 40)
 
     Lx, Ly, Lz = domain_extent.Lx, domain_extent.Ly, domain_extent.Lz
     Nx, Ny, Nz = resolution.Nx, resolution.Ny, resolution.Nz
@@ -127,11 +128,11 @@ function grid_stretching(Lz::Number, Nz::Number, refinement::Number, stretching:
 
 end
 """
-    function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number, stop_time::Number,
-                                  save_schedule::Number;  cfl = 0.75, diffusive_cfl = 0.75,
-                                  max_change = 1.2, max_Δt = 1e-1)
-Setup a DNS from `initial_conditions` that are of type `TwoLayerInitialConditions`.
-Important non-dimensional numnbers that are part of this experiment are computed and saved
+    function TLDNS_simulation_setup(tldns::TwoLayerDNS, Δt::Number, stop_time::Number,
+                                    save_schedule::Number;  cfl = 0.75, diffusive_cfl = 0.75,
+                                    max_change = 1.2, max_Δt = 1e-1)
+Setup a `Simulation` of `tldns` from `initial_conditions` that are of type `TwoLayerInitialConditions`.
+Important non-dimensional numbers that are part of this experiment are computed and saved
 to the simulation output file.
 
 ## Function arguments:
@@ -143,7 +144,7 @@ the course of a simulation;
 - `savefile` name of the file to save the data to,
 - `save_schedule` number (representing time in seconds) at which to save model output, e.g.,
 `save_schedule = 1` saves output every second;
-- `output_writer` a `Symbol` (either `:netcdf` or `:jld2`) choosing whether to save data in
+- `save_file` a `Symbol` (either `:netcdf` or `:jld2`) choosing whether to save data in
 `NetCDF` format (`.nc`) ot `JLD2` format ('.jld2).
 
 ## Keyword arguments:
@@ -158,23 +159,23 @@ there is no `Checkpointer`.
 - `max_Δt` the maximum timestep;
 - `density_reference_gp_height` for the seawater density calculation;
 - `save_velocities` defaults to `false`, if `true` model velocities will be saved to output;
-- `overwrite_existing` whether the output overwrites a file of the same name, default is
-`false`.
+- `overwrite_saved_output` whether the output overwrites a file of the same name, default is
+`true`.
 """
-function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
-                              stop_time::Number, save_schedule::Number,
-                              output_writer::Symbol=:netcdf;
-                              output_path = SIMULATION_PATH,
-                              checkpointer_time_interval = nothing,
-                              cfl = 0.75,
-                              diffusive_cfl = 0.75,
-                              max_change = 1.2,
-                              max_Δt = 1e-1,
-                              density_reference_gp_height = 0,
-                              overwrite_existing = false,
-                              save_velocities = false)
-
-    model = dns.model
+function TLDNS_simulation_setup(tldns::TwoLayerDNS, Δt::Number,
+                                stop_time::Number, save_schedule::Number;
+                                save_file = :netcdf,
+                                output_path = SIMULATION_PATH,
+                                checkpointer_time_interval = nothing,
+                                cfl = 0.75,
+                                diffusive_cfl = 0.75,
+                                max_change = 1.2,
+                                max_Δt = 1e-1,
+                                density_reference_gp_height = 0,
+                                overwrite_saved_output = nothing,
+                                save_velocities = false)
+
+    model = tldns.model
     simulation = Simulation(model; Δt, stop_time)
 
     # time step adjustments
@@ -186,62 +187,58 @@ function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
 
     # Custom saved output
     # Potential density
-    σ = seawater_density(model, geopotential_height = 0)
+    σ = seawater_density(model, geopotential_height = density_reference_gp_height)
 
     # Inferred vertical temperature diffusivity
-    T_mean = Average(T)
-    T′ = T - Field(T_mean)
-    w′ = wᶜᶜᶜ(model)
-    w′T′ = Field(w′ * T′)
-    ∫ₐw′T′ = Integral(w′T′, dims = (1, 2))
-    T_gradient = ∂z(T)
-    ∫ₐT_gradient = Integral(T_gradient, dims = (1, 2))
+    # T_mean = Average(T)
+    # T′ = T - Field(T_mean)
+    # w′ = wᶜᶜᶜ(model)
+    # w′T′ = Field(w′ * T′)
+    # ∫ₐw′T′ = Integral(w′T′, dims = (1, 2))
+    # T_gradient = ∂z(T)
+    # ∫ₐT_gradient = Integral(T_gradient, dims = (1, 2))
 
     ϵ = KineticEnergyDissipationRate(model)
     # Volume integrated TKE dissipation
     ∫ϵ = Integral(ϵ)
     # Minimum in space Kolmogorov length scale
     η_space(model) = (model.closure.ν^3 / maximum(ϵ))^(1/4)
 
-    # Volume integrated TKE dissipation
-    ∫ϵ = Integral(ϵ)
-
     # Dimensions and attributes for custom saved output
     dims = Dict("η_space" => ())
     oa = Dict(
         "σ" => Dict("longname" => "Seawater potential density calculated using TEOS-10 at $(density_reference_gp_height)dbar",
                     "units" => "kgm⁻³"),
         "η_space" => Dict("longname" => "Minimum (in space) Kolmogorov length"),
-        "∫ₐw′T′" => Dict("longname" => "Horizontally integrated vertical temperature flux (w′T′)"),
-        "∫ₐT_gradient" => Dict("longname" => "Horizontally integrated vertical temperature gradient (∂T/∂z)"),
+        # "∫ₐw′T′" => Dict("longname" => "Horizontally integrated vertical temperature flux (w′T′)"),
+        # "∫ₐT_gradient" => Dict("longname" => "Horizontally integrated vertical temperature gradient (∂T/∂z)"),
         "∫ϵ" => Dict("longname" => "Volume integrated turbulent kintetic energy dissipation")
         )
 
     # outputs to be saved during the simulation
-    outputs = Dict("S" => S, "T" => T, "η_space" => η_space, "σ" => σ, "∫ϵ" => ∫ϵ,
-                   "∫ₐw′T′" => ∫ₐw′T′, "∫ₐT_gradient" => ∫ₐT_gradient)
+    outputs = Dict("S" => S, "T" => T, "η_space" => η_space, "σ" => σ, "∫ϵ" => ∫ϵ)#=,
+                   "∫ₐw′T′" => ∫ₐw′T′, "∫ₐT_gradient" => ∫ₐT_gradient)=#
     if save_velocities
         u, v, w = model.velocities
         velocities = Dict("u" => u, "v" => v, "w" => w)
         merge!(outputs, velocities)
     end
 
-    filename = form_filename(dns, stop_time, output_writer, output_path)
-    simulation.output_writers[:outputs] = output_writer == :netcdf ?
-                                            NetCDFOutputWriter(model, outputs;
-                                                            filename,
-                                                            overwrite_existing,
+    filename = form_filename(tldns, stop_time, save_file, output_path)
+    simulation.output_writers[:outputs] = save_file == :netcdf ?
+                                            NetCDFOutputWriter(model, outputs; filename,
+                                                            overwrite_existing = overwrite_saved_output,
                                                             schedule = TimeInterval(save_schedule),
                                                             dimensions = dims,
                                                             output_attributes = oa
                                                             ) :
                                             JLD2OutputWriter(model, outputs;
                                                             filename,
                                                             schedule = TimeInterval(save_schedule),
-                                                            overwrite_existing)
+                                                            overwrite_existing = overwrite_saved_output)
 
-    non_dimensional_numbers!(simulation, dns)
-    predicted_maximum_density!(simulation, dns)
+    non_dimensional_numbers!(simulation, tldns)
+    predicted_maximum_density!(simulation, tldns)
 
     # progress reporting
     simulation.callbacks[:progress] = Callback(simulation_progress, IterationInterval(100))
@@ -252,26 +249,26 @@ function DNS_simulation_setup(dns::TwoLayerDNS, Δt::Number,
 
 end
 """
-    function form_filename(dns::TwoLayerDNS, stop_time::Number, output_writer::Symbol, output_path::AbstractString)
+    function form_filename(tldns::TwoLayerDNS, stop_time::Number, save_file, output_path)
 Create a filename for saved output based on the `profile_function`,`initial_conditions`,
 `tracer_perturbation` and length of the simulation.
 """
-function form_filename(dns::TwoLayerDNS, stop_time::Number, output_writer::Symbol, output_path::AbstractString)
+function form_filename(tldns::TwoLayerDNS, stop_time::Number, save_file, output_path)
 
-    pf_string = lowercase(string(typeof(dns.profile_function))[1:findfirst('{', string(typeof(dns.profile_function))) - 1])
-    ic_type = typeof(dns.initial_conditions)
+    pf_string = lowercase(string(typeof(tldns.profile_function))[1:findfirst('{', string(typeof(tldns.profile_function))) - 1])
+    ic_type = typeof(tldns.initial_conditions)
     ic_string = ic_type <: StableTwoLayerInitialConditions ? "stable" :
                             ic_type <: CabbelingTwoLayerInitialConditions ?
                                 "cabbeling" : ic_type <: UnstableTwoLayerInitialConditions ?
                                               "unstable" : ic_type <: IsohalineTwoLayerInitialConditions ?
                                                             "isohaline" : "isothermal"
 
-    tp_string = lowercase(string(typeof(dns.tracer_perturbation)))
+    tp_string = lowercase(string(typeof(tldns.tracer_perturbation)))
     tp_find = isnothing(findfirst('{', tp_string)) ? length(tp_string) :
                                                      findfirst('{', tp_string) - 1
     stop_time_min = stop_time / 60 ≥ 1 ? string(round(Int, stop_time / 60)) :
                                          string(round(stop_time / 60; digits = 2))
-    filetype = output_writer == :netcdf ? ".nc" : ".jld2"
+    filetype = save_file == :netcdf ? ".nc" : ".jld2"
     savefile = ic_string *"_"* pf_string *"_"* tp_string[1:tp_find] *"_"* stop_time_min * "min" * filetype
 
     # make a simulation directory if one is not present
@@ -284,11 +281,11 @@ function form_filename(dns::TwoLayerDNS, stop_time::Number, output_writer::Symbo
 
 end
 """
-    function checkpointer_setup!(simulation, model, checkpointer_time_interval)
+    function checkpointer_setup!(simulation, model, output_path, checkpointer_time_interval)
 Setup a `Checkpointer` at `checkpointer_time_interval` for a `simulation`
 """
-function checkpointer_setup!(simulation::Simulation, model::Oceananigans.AbstractModel,
-                             output_path::AbstractString, checkpointer_time_interval::Number)
+function checkpointer_setup!(simulation, model, output_path,
+                             checkpointer_time_interval::Number)
 
     dir = output_path == SIMULATION_PATH ? CHECKPOINT_PATH :
                                            joinpath(output_path, "model_checkpoints/")

test/initialconditions_test.jl

@@ -1,8 +1,8 @@
 # Need to write these to check things are set at the right level in the right field
 diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
 resolution = (Nx = 10, Ny = 10, Nz = 500)
-model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
-            reference_density = REFERENCE_DENSITY)
+model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+                reference_density = REFERENCE_DENSITY)
 
 ## set initial conditions, currently there are four options available in this submodule
 initial_conditions = StableTwoLayerInitialConditions(0, 0, 0, 0, 0, 0)

test/kernelfunction_test.jl

@@ -7,8 +7,8 @@ using GibbsSeaWater
 diffusivities = (ν = 1e-4, κ = (S = 1e-5, T = 1e-5))
 resolution = (Nx = 10, Ny = 10, Nz = 100)
 
-model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
-            reference_density = REFERENCE_DENSITY)
+model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+                reference_density = REFERENCE_DENSITY)
 
 z = znodes(model.grid, Center())
 

test/runtests.jl

@@ -63,7 +63,7 @@ include("initialconditions_test.jl")
 
     for tb ∈ tracer_profile_perturbations
 
-        model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+        model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
                     reference_density = REFERENCE_DENSITY)
         dns = TwoLayerDNS(model, profile_function, initial_conditions, tracer_perturbation = tb)
         set_two_layer_initial_conditions!(dns)
@@ -80,7 +80,7 @@ end
 
     for tb ∈ tracer_blob_perturbations
 
-        model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+        model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
                     reference_density = REFERENCE_DENSITY)
         dns = TwoLayerDNS(model, profile_function, initial_conditions, tracer_perturbation = tb)
         set_two_layer_initial_conditions!(dns)
@@ -94,7 +94,7 @@ end
 
     for tn ∈ tracer_noise_perturbations
 
-        model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+        model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
                     reference_density = REFERENCE_DENSITY)
         dns = TwoLayerDNS(model, profile_function, initial_conditions, initial_noise = tn)
         set_two_layer_initial_conditions!(dns)
@@ -104,7 +104,7 @@ end
     end
     for tnv ∈ tracer_noise_perturbations_vec
 
-        model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+        model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
                     reference_density = REFERENCE_DENSITY)
         dns = TwoLayerDNS(model, profile_function, initial_conditions, initial_noise = tnv)
         set_two_layer_initial_conditions!(dns)
@@ -117,7 +117,7 @@ end
 
 @testset "Step change" begin
 
-    model = DNS(architecture, DOMAIN_EXTENT, resolution, diffusivities;
+    model = DNSModel(architecture, DOMAIN_EXTENT, resolution, diffusivities;
                 reference_density = REFERENCE_DENSITY)
     profile_function = StepChange(z[depth_idx])
     initial_conditions = TwoLayerInitialConditions(34.551, -1.5, 34.7, 0.5)