Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow users to specify the domain when creating a `RegularCartes… #464

Merged
merged 15 commits into from
Oct 23, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion benchmark/benchmark_static_ocean.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ for arch in archs, float_type in float_types, N in Ns
Nx, Ny, Nz = N
Lx, Ly, Lz = 1, 1, 1

model = Model(architecture=arch, float_type=float_type, grid=RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)))
model = Model(architecture=arch, float_type=float_type, grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)))
time_step!(model, Ni, 1)

bname = benchmark_name(N, "", arch, float_type)
Expand Down
5 changes: 2 additions & 3 deletions benchmark/benchmark_tracers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,10 +57,10 @@ for arch in archs, test_case in test_cases

bname = benchmark_name(N, "$na active + $(lpad(np, 2)) passive", arch, FT)
@printf("Running benchmark: %s...\n", bname)

model = Model(architecture = arch,
float_type = FT,
grid = RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
buoyancy = na2buoyancy(na),
tracers = tracers)
time_step!(model, Ni, 1)
Expand All @@ -72,4 +72,3 @@ end

print_timer(timer, title="Tracer benchmarks", sortby=:name)
println("")

6 changes: 3 additions & 3 deletions examples/internal_wave.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,17 +56,17 @@ a(x, z) = A * exp( -( (x - x₀)^2 + (z - z₀)^2 ) / 2δ^2 )
u₀(x, y, z) = a(x, z) * U * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * V * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * W * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z
b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z

# We are now ready to instantiate our model on a uniform grid.
# We give the model a constant rotation rate with background vorticity `f`,
# use temperature as a buoyancy tracer, and use a small constant viscosity
# and diffusivity to stabilize the model.

model = Model(
grid = RegularCartesianGrid(N=(Nx, 1, Nx), L=(Lx, Lx, Lx)),
grid = RegularCartesianGrid(size=(Nx, 1, Nx), length=(Lx, Lx, Lx)),
closure = ConstantIsotropicDiffusivity(ν=1e-6, κ=1e-6),
coriolis = FPlane(f=f),
coriolis = FPlane(f=f),
tracers = :b,
buoyancy = BuoyancyTracer()
)
Expand Down
2 changes: 1 addition & 1 deletion examples/netcdf_ouput_example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Nx, Ny, Nz = 16, 16, 16 # No. of grid points in x, y, and z, respectively.
Lx, Ly, Lz = 100, 100, 100 # Length of the domain in x, y, and z, respectively (m).
tf = 5000 # Length of the simulation (s)

model = Model(grid=RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
model = Model(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
closure=ConstantIsotropicDiffusivity())

# Add a cube-shaped warm temperature anomaly that takes up the middle 50%
Expand Down
14 changes: 7 additions & 7 deletions examples/ocean_convection_with_plankton.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,20 +6,20 @@
# * to set boundary conditions;
# * to defined and insert a user-defined forcing function into a simulation.
# * to use the `TimeStepWizard` to manage and adapt the simulation time-step.
#
#
# To begin, we load Oceananigans, a plotting package, and a few miscellaneous useful packages.

using Oceananigans, PyPlot, Random, Printf

# ## Parameters
#
# We choose a modest two-dimensional resolution of 128² in a 64² m² domain ,
# implying a resolution of 0.5 m. Our fluid is initially stratified with
# implying a resolution of 0.5 m. Our fluid is initially stratified with
# a squared buoyancy frequency
#
# $$ N\^2 = 10^{-5} \, \mathrm{s^{-2}} $$
#
# and a surface buoyancy flux
# and a surface buoyancy flux
#
# $$ Q_b2 = 10^{-8} \, \mathrm{m^3 \, s^{-2}} $$
#
Expand All @@ -39,7 +39,7 @@ end_time = 1day
# Create boundary conditions. Note that temperature is buoyancy in our problem.
#

buoyancy_bcs = HorizontallyPeriodicBCs( top = BoundaryCondition(Flux, Qb),
buoyancy_bcs = HorizontallyPeriodicBCs( top = BoundaryCondition(Flux, Qb),
bottom = BoundaryCondition(Gradient, N²))

# ## Define a forcing function
Expand All @@ -51,10 +51,10 @@ buoyancy_bcs = HorizontallyPeriodicBCs( top = BoundaryCondition(Flux, Qb),
growth_and_decay = SimpleForcing((x, y, z, t) -> exp(z/16) - 1)

## Instantiate the model
model = Model(
grid = RegularCartesianGrid(N = (Nz, 1, Nz), L = (Lz, Lz, Lz)),
model = Model(
grid = RegularCartesianGrid(size = (Nz, 1, Nz), length = (Lz, Lz, Lz)),
closure = ConstantIsotropicDiffusivity(ν=1e-4, κ=1e-4),
coriolis = FPlane(f=1e-4),
coriolis = FPlane(f=1e-4),
tracers = (:b, :plankton),
buoyancy = BuoyancyTracer(),
forcing = ModelForcing(plankton=growth_and_decay),
Expand Down
54 changes: 27 additions & 27 deletions examples/ocean_wind_mixing_and_convection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,30 +18,30 @@ using Oceananigans, PyPlot, Random, Printf
# Here we use an isotropic, cubic grid with `Nz` grid points and grid spacing
# `Δz = 1` meter. We specify fluxes of heat, momentum, and salinity via
#
# 1. A temperature flux `Qᵀ` at the top of the domain, which is related to heat flux
# by `Qᵀ = Qʰ / (ρ₀ * cᴾ)`, where `Qʰ` is the heat flux, `ρ₀` is a reference density,
# and `cᴾ` is the heat capacity of seawater. With a reference density
# `ρ₀ = 1026 kg m⁻³`and heat capacity `cᴾ = 3991`, our chosen temperature flux of
# `Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹` corresponds to a heat flux of `Qʰ = 204.7 W m⁻²`, a
# 1. A temperature flux `Qᵀ` at the top of the domain, which is related to heat flux
# by `Qᵀ = Qʰ / (ρ₀ * cᴾ)`, where `Qʰ` is the heat flux, `ρ₀` is a reference density,
# and `cᴾ` is the heat capacity of seawater. With a reference density
# `ρ₀ = 1026 kg m⁻³`and heat capacity `cᴾ = 3991`, our chosen temperature flux of
# `Qᵀ = 5 × 10⁻⁵ K m⁻¹ s⁻¹` corresponds to a heat flux of `Qʰ = 204.7 W m⁻²`, a
# relatively powerful cooling rate.
#
# 2. A velocity flux `Qᵘ` at the top of the domain, which is related
# to the `x` momentum flux `τˣ` via `τˣ = ρ₀ * Qᵘ`, where `ρ₀` is a reference density.
# Our chosen value of `Qᵘ = -2 × 10⁻⁵ m² s⁻²` roughly corresponds to atmospheric winds
# of `uᵃ = 2.9 m s⁻¹` in the positive `x`-direction, using the parameterization
# of `uᵃ = 2.9 m s⁻¹` in the positive `x`-direction, using the parameterization
# `τ = 0.0025 * |uᵃ| * uᵃ`.
#
# 3. An evaporation rate `evaporation = 10⁻⁷ m s⁻¹`, or approximately 0.1 millimeter per
# 3. An evaporation rate `evaporation = 10⁻⁷ m s⁻¹`, or approximately 0.1 millimeter per
# hour.
#
# Finally, we use an initial temperature gradient of `∂T/∂z = 0.005 K m⁻¹`,
# which implies an iniital buoyancy frequency `N² = α * g * ∂T/∂z = 9.8 × 10⁻⁶ s⁻²`
# with a thermal expansion coefficient `α = 2 × 10⁻⁴ K⁻¹` and gravitational acceleration
# `g = 9.81 s⁻²`. Note that, by default, the `SeawaterBuoyancy` model uses a gravitational
# with a thermal expansion coefficient `α = 2 × 10⁻⁴ K⁻¹` and gravitational acceleration
# `g = 9.81 s⁻²`. Note that, by default, the `SeawaterBuoyancy` model uses a gravitational
# acceleration `gᴱᵃʳᵗʰ = 9.80665 s⁻²`.

Nz = 48 # Number of grid points in x, y, z
Δz = 1.0 # Grid spacing in x, y, z (meters)
Δz = 1.0 # Grid spacing in x, y, z (meters)
Qᵀ = 5e-5 # Temperature flux at surface
Qᵘ = -2e-5 # Velocity flux at surface
∂T∂z = 0.005 # Initial vertical temperature gradient
Expand Down Expand Up @@ -70,17 +70,17 @@ S_bcs = HorizontallyPeriodicBCs(top = BoundaryCondition(Flux, Qˢ))

# ## Model instantiation
#
# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`,
# using a `FPlane` model for rotation (constant rotation rate), a linear equation
# of state for temperature and salinity, the Anisotropic Minimum Dissipation closure
# We instantiate a horizontally-periodic `Model` on the CPU with on a `RegularCartesianGrid`,
# using a `FPlane` model for rotation (constant rotation rate), a linear equation
# of state for temperature and salinity, the Anisotropic Minimum Dissipation closure
# to model the effects of unresolved turbulence, and the previously defined boundary
# conditions for `u`, `T`, and `S`. We also pass the evaporation rate to the container
# model.parameters for use in the boundary condition function that calculates the salinity
# flux.

model = Model(
architecture = CPU(),
grid = RegularCartesianGrid(N=(Nz, Nz, Nz), L=(Δz*Nz, Δz*Nz, Δz*Nz)),
grid = RegularCartesianGrid(size=(Nz, Nz, Nz), length=(Δz*Nz, Δz*Nz, Δz*Nz)),
coriolis = FPlane(f=f),
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
closure = AnisotropicMinimumDissipation(),
Expand All @@ -95,14 +95,14 @@ model = Model(
# ```
#
# To change the `architecture` to `GPU`, replace the `architecture` keyword argument with
#
#
# ```julia
# architecture = GPU()
# ```
#
# ## Initial conditions
#
# Out initial condition for temperature consists of a linear stratification superposed with
# Out initial condition for temperature consists of a linear stratification superposed with
# random noise damped at the walls, while our initial condition for velocity consists
# only of random noise.

Expand All @@ -120,24 +120,24 @@ set!(model, u=u₀, w=u₀, T=T₀, S=35)
# ## Set up output
#
# We set up an output writer that saves all velocity fields, tracer fields, and the subgrid
# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument
# to `JLD2OutputWriter` indicates that output will be saved in
# turbulent diffusivity associated with `model.closure`. The `prefix` keyword argument
# to `JLD2OutputWriter` indicates that output will be saved in
# `ocean_wind_mixing_and_convection.jld2`.

## Create a NamedTuple containing all the fields to be outputted.
fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))

## Instantiate a JLD2OutputWriter to write fields.
field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
field_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output); interval=hour/4,
prefix="ocean_wind_mixing_and_convection", force=true)

## Add the output writer to the models `output_writers`.
model.output_writers[:fields] = field_writer

# ## Running the simulation
#
# To run the simulation, we instantiate a `TimeStepWizard` to ensure stable time-stepping
# with a Courant-Freidrichs-Lewy (CFL) number of 0.2.
# with a Courant-Freidrichs-Lewy (CFL) number of 0.2.

wizard = TimeStepWizard(cfl=0.2, Δt=1.0, max_change=1.1, max_Δt=5.0)

Expand All @@ -153,11 +153,11 @@ fig, axs = subplots(ncols=3, figsize=(12, 5))
"""
makeplot!(axs, model)

Make a triptych of x-z slices of vertical velocity, temperature, and salinity
Make a triptych of x-z slices of vertical velocity, temperature, and salinity
associated with `model` in `axs`.
"""
function makeplot!(axs, model)
jhalf = floor(Int, model.grid.Nz/2)
jhalf = floor(Int, model.grid.Nz/2)

## Coordinate arrays for plotting
xC = repeat(model.grid.xC, 1, model.grid.Nz)
Expand All @@ -169,7 +169,7 @@ function makeplot!(axs, model)
pcolormesh(xC, zF, data(model.velocities.w)[:, jhalf, :])
xlabel("\$ x \$ (m)"); ylabel("\$ z \$ (m)")

sca(axs[2]); cla()
sca(axs[2]); cla()
title("Temperature")
pcolormesh(xC, zC, data(model.tracers.T)[:, jhalf, :])
xlabel("\$ x \$ (m)")
Expand All @@ -196,9 +196,9 @@ while model.clock.time < end_time
walltime = @elapsed time_step!(model, 10, wizard.Δt)

## Print a progress message
@printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
wmax(model), prettytime(walltime))
@printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
wmax(model), prettytime(walltime))

model.architecture == CPU() && makeplot!(axs, model)
end
2 changes: 1 addition & 1 deletion examples/simple_diffusion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ using PyPlot, Printf
# `Model` constructor:

model = Model(
grid = RegularCartesianGrid(N = (1, 1, 128), L = (1, 1, 1)),
grid = RegularCartesianGrid(size = (1, 1, 128), length = (1, 1, 1)),
closure = ConstantIsotropicDiffusivity(κ = 1.0)
)

Expand Down
4 changes: 2 additions & 2 deletions examples/two_dimensional_turbulence.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ using Oceananigans.TurbulenceClosures: ∂x_faa, ∂y_afa
# We instantiate the model with a simple isotropic diffusivity

model = Model(
grid = RegularCartesianGrid(N=(128, 128, 1), L=(2π, 2π, 2π)),
buoyancy = nothing,
grid = RegularCartesianGrid(size=(128, 128, 1), length=(2π, 2π, 2π)),
buoyancy = nothing,
tracers = nothing,
closure = ConstantIsotropicDiffusivity(ν=1e-3, κ=1e-3)
)
Expand Down
2 changes: 1 addition & 1 deletion src/fields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ a function with arguments `(x, y, z)`, or any data type for which a
Example
=======
```julia
model = Model(grid=RegularCartesianGrid(N=(32, 32, 32), L=(1, 1, 1))
model = Model(grid=RegularCartesianGrid(size=(32, 32, 32), length=(1, 1, 1))

# Set u to a parabolic function of z, v to random numbers damped
# at top and bottom, and T to some silly array of half zeros,
Expand Down
Loading