Skip to content

Commit

Permalink
Allow users to specify the domain when creating a `RegularCartes… (Cl…
Browse files Browse the repository at this point in the history
…iMA#464)

Allow users to specify the domain when creating a `RegularCartesianGrid`
  • Loading branch information
ali-ramadhan authored Oct 23, 2019
2 parents e98cb43 + cab311e commit 3f3ef52
Show file tree
Hide file tree
Showing 24 changed files with 197 additions and 170 deletions.
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

0 comments on commit 3f3ef52

Please sign in to comment.